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A  methodology  has  been  developed  to  simulate  fluid  flow,  heat  transfer,  free  surface, 
and  phase  change  dynamics  in  the  context  of  solidification  processing.  Special  attention  has 
been  given  to  float  zone  processing  of  single  crystals  due  to  its  increasing  technological 
importance  in  growing  high  quality  single  crystals  and  due  to  the  complex  physical 
phenomena  involved  -  free  surface  flow,  thermocapillary  and  buoyancy  driven  convection, 
phase  change  and  the  nonlinear  coupling  amongst  these  physical  phenomena.  The  SIMPLE 
algorithm  in  body  fitted  coordinates  has  been  chosen  for  fluid  flow  and  heat  transfer 
computations.  Issues  of  accuracy  and  implementation  of  convection  schemes  have  been 
addressed  in  the  context  of  buoyancy  driven  convection.  Accurate  simulations  of  transient 
natural  convection  around  an  enclosed  vertical  channel  have  been  conducted  since  it  is  an 
important  heat  transfer  mechanism  in  solidification  as  well  as  a  host  of  other  applications 
such  as  electronics  cooling.  The  issues  of  existence,  multiplicity  and  non-uniqueness  of  float 
zone  menisci  have  been  dealt  with  in  the  context  of  a  free  energy  minimization  principle. 
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Steady,  thermocapillary  convection  has  been  simulated  in  a  liquid  bridge  which  yields  useful 
information  on  heat  transfer  characteristics  in  a  float  zone  configuration.  Both  single  phase 
flows  and  flows  with  phase  change  have  been  computed  in  this  study.  Alternative 
formulations  of  enthalpy  porosity  models  have  been  developed  and  assessed  for  the 
computation  of  phase  change  dynamics  in  the  presence  of  melt  convection.  These  methods 
have  been  applied  to  compute  and  assess  the  impact  of  interacting  natural  and 
thermocapillary  convection  on  solid-liquid  interfaces  during  float  zone  growth.  A  moving 
boundary  algorithm  has  been  developed  to  simulate  the  nonlinear  coupling  between  the 
unknown  shape  and  location  of  the  free  surface,  the  unknown  shape  and  location  of  the 
solid-liquid  interfaces,  and  the  convective  heat  transfer  processes  in  the  float  zone.  This 
computational  model  has  also  been  applied  to  assess  the  impact  of  both,  normal  gravity  as 
well  as  microgravity  conditions,  on  the  dynamics  within  the  float  zone.  It  has  been  found  that 
whereas  buoyancy  induced  convection  is  negligible  under  microgravity  conditions, 
significant  levels  of  thermocapillary  convection  exist  and  have  a  profound  impact  on  the 
solid-liquid  interface  morphology  at  the  macroscopic  scales.  It  has  also  been  found  that 
under  normal  gravity  conditions,  the  counteracting  effect  of  natural  and  thermocapillary 
convection  can  actually  result  in  reduced  convection  strength  in  the  lower  half  of  the  zone 
and  lead  to  a  flatter  solid-liquid  interface  shape.  The  calculation  procedures  developed  in 
this  study  constitute  a  predictive  capability  for  float  zone  single  crystal  growth. 
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CHAPTER  1 
INTRODUCTION 


1 . 1  Motivation 

The  goal  of  this  research  is  to  develop  and  employ  computational  tools  to  model  phase 
change  and  convective  heat  transfer  in  materials  processing  applications.  Special  emphasis 
is  placed  on  problems  dealing  with  solidification  and  crystal  growth.  Typical  applications 
include  growth  of  optical  and  semiconductor  crystals  from  the  melt  as  well  as  intermetallics 
for  high  temperature  structural  applications.  Solidification  phenomena  occurring  during 
materials  processing  are  dominated  by  convective  heat  and  mass  transport  in  the  melt,  and 
more  importantly,  are  characterized  by  the  presence  of  one  or  more  moving  interfaces 
constituting  moving  boundaries.  These  interfaces  are  discontinuities  separating  different 
flow  regions;  they  may  be  phase  change  interfaces  such  as  those  between  the  melt  and  the 
solid  or  they  may  be  free  surfaces  separating  the  melt  and  the  ambient  liquid  or  gas.  The 
nonlinear  coupling  between  these  different  phenomena  is  the  major  source  of  modeling  and 
computational  difficulties. 

Modeling  of  fluid  flow  and  heat  transfer  is  a  matter  of  particular  interest  in  this  study. 
These  transport  processes  control  the  quality  of  the  grown  crystal  or  solidified  product 
(Brown  1988,  Huppert  1990,  Davis  1990).  For  example,  high  thermal  gradients  in  the 
growing  crystal  can  lead  to  large  thermoelastic  stresses  that  can  cause  cracking  or  fracture 
(Viskanta  1990).  Solute  rejection  at  the  growing  interface  can  lead  to  constitutional 
undercooling  giving  rise  to  morphological  instabilities  of  the  interface  (Mullins  and  Sekerka 
1964,  Glicksman  et  al.  1986). 
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1.2  Computational  Methodology 

Fluid  convection  can  arise  in  solidification  systems  due  to  a  variety  of  factors.  Density 
differences  in  the  melt  caused  by  thermal  and  compositional  inhomogeneities  give  rise  to 
buoyancy  induced  convection  (Langlois  1985).  Marangoni  convection  arises  out  of  shear 
stress  induced  by  thermal  and  composition  gradients  along  free  surfaces.  Density  differences 
between  the  melt  and  the  solid  can  also  be  an  important  source  of  convection  depending  on 
the  growth  configuration.  In  addition,  forced  convection  can  be  present  due  to  crystal  and/or 
crucible  rotation  which  is  usually  done  in  order  to  even  out  azimuthal  inhomogeneities  in 
the  melt  temperature  and  composition  (Langlois  1985,  Hurle  1993). 

In  this  work,  attention  is  focused  on  thermally  induced  buoyancy  driven  and 
thermocapillary  convection,  both,  with  and  without  the  presence  of  phase  change.  Whereas 
a  major  interest  is  in  problems  dealing  with  crystal  growth,  natural  and  thermocapillary 
convection  without  phase  change  have  been  considered  as  a  natural  starting  point  to  develop 
and  assemble  the  necessary  computational  capabilities.  The  SIMPLE  algorithm  for 
incompressible  flow  computation  (Patankar  1980)  extended  to  treat  complex  geometry  by 
using  a  body  fitted  coordinate  system  (Shyy  1 994)  is  chosen  as  the  underlying  computational 
technique.  It  is  a  robust  and  proven  method,  and  the  sequential  solution  technique  built  into 
the  method  enables  additional  physics  to  be  incorporated  without  the  need  for  substantial 
modification  to  the  algorithm.  Thus,  the  computational  technique  can  be  extended  in  stages 
that  include  more  complex  physical  phenomena  which  increases  the  level  of  confidence  in 
the  computational  technique. 

As  a  starting  point  and  to  develop  the  necessary  numerical  techniques,  problems 
involving  natural  and  thermocapillary  convection,  but  without  phase  change,  have  been 
investigated.  The  SIMPLE  algorithm  has  been  implemented  and  tested  against  an  available 
benchmark  solution  for  natural  convection  of  air  in  a  square  cavity.  In  chapter  3,  this  problem 
is  also  used  to  characterize  the  relative  accuracy  of  various  convection  schemes  such  as  the 
second  order  central  difference  scheme,  the  second  order  upwind  scheme  etc.  for  the 
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momentum  and  energy  equations.  The  computational  technique  has  also  been  extended  to 
handle  the  unsteady  terms  in  the  Navier-Stokes  equations.  The  computational  technique  is 
then  applied  to  model  transient  natural  convection  and  heat  transfer  in  a  configuration  used 
for  electronics  cooling.  Chapter  3  also  describes  the  extension  of  the  computational 
procedure  to  simulate  steady  thermocapillary  convection  in  a  liquid  bridge  of  molten  NiAl 
-  an  intermetallic  that  is  being  considered  for  high  temperature  structural  application  in 
future  generation  of  aircraft  engines  and  other  aerospace  applications.  The  model  is  based 
on  a  float  zone  configuration  currently  in  use  for  growing  single  crystals  of  NiAl. 

In  chapter  4,  the  computational  technique  is  further  extended  to  model  phase  change  and 
solidification  processes.  The  enthalpy  method  (Crank  1984,  Shyy  1994,  Shyy  et.  al.  1995) 
is  particularly  suitable  for  modeling  phase  change  processes  at  macroscopic  scales  and  is 
chosen  for  conducting  the  solidification  simulations.  The  advantage  of  the  enthalpy  method 
is  that  it  is  a  fixed  grid  technique  which  obviates  the  need  for  a  moving  grid  and  explicit 
interface  tracking.  However,  the  volume  averaging  procedure  that  is  implicit  in  this 
technique  smears  out  the  interface  information  over  a  thin  mushy  zone  implying  that 
curvature  dependent  phenomena  and  interface  microstructure  cannot  be  captured.  A 
one-dimensional  heat  conduction  problem  (Stefan  problem)  is  chosen  as  the  test  bed  to 
develop  a  particularly  efficient  implementation  of  the  enthalpy  method  and  to  characterize 
its  accuracy  and  performance.  The  enthalpy  method  is  also  applied  to  a  two  dimensional 
melting  problem  dominated  by  natural  convection  heat  transfer  and  highly  nonlinear 
interfacial  dynamics.  Qualitative  agreement  of  the  computed  interface  shapes  with 
experimental  data  is  established. 

The  enthalpy  method  is  then  applied  to  simulate  practical  crystal  growth  configurations. 
The  float  zone  model  considered  in  chapter  3  is  extended,  in  chapter  4,  to  take  into  account 
the  nonlinear  coupling  of  interacting  natural  and  thermocapillary  convection  with  the 
melt/crystal  and  melt/feed  interfaces.  However,  the  free  surface  is  held  stationary  in  this 
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study  in  order  to  clearly  isolate  the  interaction  of  convection  heat  transfer  and  phase  change 
dynamics.  This  restriction  is  removed  in  the  studies  taken  up  in  Chapter  5. 


The  rest  of  this  work  deals  with  the  extension  of  the  computational  techniques 
developed,  to  compute  the  unknown  shape  and  location  of  the  free  surface  coupled  with 
phase  change  and  convective  heat  transfer.  The  float  zone  growth  technique  is  taken  as  the 
prototype  application.  The  float  zone  process  consists  of  the  growing  crystal  and  the  feed 
separated  by  a  liquid  bridge  of  the  molten  material  as  shown  in  Fig.  1.1.  An  adequate 
simulation  involves  consideration  of  the  nonlinear  coupling  of  the  melt/crystal  and 
feed/crystal  interfaces,  the  shape  and  motion  of  the  free  surface  (meniscus)  between  the  melt 
and  ambient  gas  or  encapsulant,  thermocapillary  and  buoyancy  driven  convection  in  the 
melt.  In  order  to  build  an  effective  simulation  capability,  a  moving  grid  technique  similar 
to  the  one  developed  in  Shyy  et.  al.  (1995)  is  developed  to  track  the  motion  of  the  free 
surface.  The  shape  and  location  of  the  free  surface  are  determined  by  integrating  the  dynamic 


1.3  Float  Zone  Crystal  Growth 
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Figure  1.1.  Schematic  of  the  float  zone  process.  The  melt  zone  is  sustained 
by  the  surface  tension  between  the  melt  and  the  encapsulant.  The  heating  profile 
is  schematically  illustrated  in  the  right. 
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form  of  the  Young-Laplace  equation.  The  Young-Laplace  equation  describes  the  balance 
between  the  surface  tension  force  and  the  pressure  difference  across  the  free  surface. 

The  float  zone  technique  is  one  of  the  candidate  methods  for  growing  single  crystals. 
The  main  advantage  of  the  float  zone  process  is  that  it  is  a  containerless  process  and  thus  one 
main  source  of  contamination  is  removed.  Another  advantage  is  the  reduction  of  thermal 
stresses  which  are  caused  by  differential  thermal  expansion  between  the  crystal  and  the 
crucible.  On  earth,  it  is  mainly  limited  by  the  size  of  the  melt  zone  that  can  be  achieved, 
since  the  hydrostatic  pressure  of  the  melt  zone  is  balanced  by  the  surface  tension  between 
the  melt  and  the  ambient  fluid.  Another  problem  is  the  loss  of  volatile  components  from  the 
melt,  which  can  be  circumvented  by  the  use  of  encapsulating  fluids.  Liquid  encapsulants 
can  also  increase  the  maximum  size  of  the  melt  zone  that  can  be  achieved.  The  formation 
of  the  meniscus  during  the  melting  process  is  a  major  concern  in  the  successful  operation  of 
the  float  zone  technique  and  this  aspect  will  be  examined  in  detail  in  chapter  2. 

Float  zone  simulations  have  been  carried  out  by  many  researchers  at  different  levels  of 
complexity.  Kobayashi  (1984)  computed  thermocapillary  convection  under  zero  gravity 
conditions  for  two  different  Prandtl  numbers,  with  and  without  crystal  rotation  while 
neglecting  the  deformation  of  the  free  surface.  Kobayashi  (1988)  extended  the  study  to 
include  buoyancy  induced  convection.  Young  (1990)  has  presented  asymptotic  solutions  for 
the  temperature,  concentration,  and  interface  shapes  in  the  limit  of  small  aspect  ratio  and 
weak  surface  tension  for  a  two  dimensional  float  zone  established  in  a  vertical  sheet.  Neitzel 
et.  al.  ( 1 990)  have  computed  thermocapillary  convection  in  a  float  zone  with  a  computed  free 
surface  shape.  Chen  et.  al.  (1994)  have  carried  out  numerical  simulation  for  low  and  high 
Prandtl  number  fluids  to  study  melt  convection  and  the  power  input  required  to  form  a 
molten  zone.  However,  they  have  constrained  the  free  surface  to  a  cylindrical  shape  and  have 
not  considered  the  effect  of  free  surface  deformation  and  curvature.  Chu  and  Chen  (1995) 
have  obtained  numerical  solutions  for  float  zone  growth  of  Molybdenum  (Mo)  crystals  with 
an  imposed  heat  input  profile  rather  than  an  imposed  temperature  profile.  Again,  they 
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imposed  the  constraint  of  a  fixed  cylindrical  free  surface  and  concentrated  their  efforts  on 
the  details  of  the  flow  structure.  Duranceau  and  Brown  (1986)  developed  a  steady-state, 
thermal-capillary  model  where  the  free  surface  was  computed  as  part  of  the  solution. 
However,  their  model  considered  only  heat  conduction  but  not  melt  convection.  Lan  and 
Kou  (1990a,  1990b)  have  carried  out  coupled  calculations  of  steady  thermocapillary  flow 
and  natural  convection  under  normal  and  microgravity  conditions  for  an  NaNC»3  float  zone. 
However,  they  too  confined  their  calculations  to  an  imposed  cylindrical  shape  of  the  free 
surface. 

Duranceau  and  Brown  (1986)  developed  a  steady-state,  thermal-capillary  model  where 
the  free  surface  was  computed  as  part  of  the  solution.  However,  their  model  considered  only 
heat  conduction  but  not  melt  convection.  Lan  and  Kou  (1991a,  1991b,  1992a)  carried  out 
computations  of  float  zone  growth  of  NaNC>3  crystals  where  the  heater  was  in  physical 
contact  with  the  melt.  Their  intention  was  to  grow  crystals  of  uniform  cross-section  and 
improved  surface  smoothness  as  well  as  melt  zone  stability  at  the  cost  of  increased  risk  of 
melt  contamination.  This  arrangement  prevented  the  deformation  of  the  free  surface  and 
suppressed  thermocapillary  convection  to  a  large  extent.  They  also  considered  cases  where 
the  melt  was  partially  exposed  to  the  ambient  gas  thus  permitting  weak  thermocapillary 
convection  but  not  free  surface  deformation. 

Lan  and  Kou  (1991c,  1992b)  have  also  carried  out  a  series  of  studies  on  float  zones  with 
rotation  and  counterrotation  of  the  feed  and  crystal.  The  have  found  that  due  to  the  resulting 
melt  rotation  the  height  of  the  zone  can  be  reduced  significantly  by  reducing  the  power  input. 
The  forced  convection  in  the  melt  due  to  rotation  helps  to  prevent  freezing  of  the  melt  at  the 
centerline.  This  is  especially  useful  for  non-electrically  conducting  materials  such  as 
NaN03  since  electromagnetic  levitation  forces  cannot  be  used  to  stabilize  the  zone. 

Chan  and  Choi  (1992)  and  Chan  (1995)  have  modeled  the  convective  heat  transfer, 
solidification  and  free  surface  deformation  during  float  zone  growth  of  Titanium  Carbide 
(TiC)  by  induction  heating.  Muhlbauer  et.  al.  (1995)  have  simulated  the  growth  of  large 
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Silicon  crystals  by  the  needle-eye  technique  with  the  unknown  interface  shapes  determined 
as  a  coupled  thermal-electromagnetic-hydrodynamic  problem.  For  their  problem,  buoyancy 
induced  convection  and  forced  convection  due  to  crystal  rotation  dominated  thermocapillary 
convection  in  the  float  zone.  The  effect  of  dopant  segregation  in  float  zone  growth  of  KNO3 
doped  NaNC>3  crystals  has  been  considered  by  Lan  and  Kou  (1993a)  at  zero  gravity.  Also, 
in  a  related  work,  Lan  and  Kou  (1993b)  have  examined  dopant  segregation  under 
microgravity  and  with  crystal  rotation.  In  these  studies,  both  negative  and  positive 
segregation  coefficients  have  been  considered.  The  have  found  that  the  shape  of  the 
melt/crystal  interface  is  significantly  affected  by  the  radial  dopant  segregation,  especially  at 
the  centerline.  Also,  with  increasing  strength  of  thermocapillary  convection,  the  segregation 
decreases.  A  similar  result  was  obtained  with  fast  counterrotation  of  the  crystal.  Lan  (1994) 
computed  steady,  convective  heat  transfer  and  interface  shapes  in  a  float  zone  using 
Newton's  iterations.  This  study  employed  a  mass  conservation  scheme  to  determine  the  free 
surface  shape  as  opposed  to  a  fixed  contact  angle  condition. 

Transient  computations  of  g-jitter  induced  convection  in  a  Bismuth  Silicate  (BSO)  float 
zone  have  been  studied  by  Chen  et.  al.  (1994).  However,  they  constrained  the  free  surface 
to  a  nondeformable  cylindrical  shape  and  imposed  a  fixed,  parabolic  temperature 
distribution  on  the  free  surface  and  did  not  consider  phase  change/solidification  effects. 
Grugel  et.  al.  (1994)  have  carried  out  experimental  investigations  of  vibration  on  the 
solidification  characteristics  in  a  NaN03  float  zone.  They  found  that  vibration  suppressed 
thermocapillary  convection  resulting  in  flatter  interfaces  as  evidenced  by  the  improved 
uniformity  of  the  microstructural  morphology  presumably  due  to  the  opposing  effects  of 
thermocapillary  and  buoyancy  driven  convection  set  up  the  by  effect  of  body  forces  due  to 
vibration.  Hu  et.  al.  (1994)  have  carried  out  experiments  on  the  onset  of  oscillatory 
thermocapillary  flow  in  a  half  zone  to  determine  the  relationship  between  the  liquid  bridge 
volume  and  the  critical  Marangoni  number  measured  in  terms  of  the  critical  temperature 
difference.  Tang  and  Hu  (1994)  carried  out  corresponding  numerical  simulations  on  a  two 
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dimensional  half  zone.  Lee  et.  al.  (1995)  have  performed  experiments  to  investigate 
oscillatory  thermocapillary  flow.  They  found  that  three  dimensional  oscillatory  flow  appears 
when  the  Marangoni  number  exceeds  a  critical  value.  The  flow  modes  were  found  to  be 
strongly  dependent  on  the  aspect  ratio  rather  than  on  the  Marangoni  number.  The  flow 
oscillations,  thermal  oscillations  and  the  free  surface  deformations  had  the  same  frequencies. 
Yeckel  et.  al.  (1995)  have  simulated  heat  and  momentum  transport  in  the  vertical  float  zone 
refinement  of  two  dimensional  silicon  sheets.  They  found  that  thermocapillary  convection 
dominates  over  buoyancy  driven  flow  and  the  solid-liquid  interfaces  depart  considerably 
from  a  planar  shape.  Levenstam  et.  al.  (1996)  have  carried  out  experimental  and  theoretical 
studies  to  determine  the  critical  Reynolds  number  for  transition  to  oscillatory  flow  in  a 
silicon  float  zone.  They  found  that  as  the  Reynolds  number  increases,  the  flow  transitions 
from  a  steady  two  dimensional  state  to  a  three  dimensional  state  and  then  to  an  unsteady  flow. 

Thus,  many  researchers  have  worked  towards  the  goal  of  float  zone  simulations  with 
varying  degrees  of  numerical  sophistication,  physical  complexity  and  simplifications 
regarding  the  free  surface  shape  etc.  The  simulation  capability  developed  in  this  work  aims 
at  a  more  complete  level  of  modeling  capability.  A  combination  of  higher  grid  resolution, 
second  order  convection  schemes,  strict  convergence  criteria,  and  moving  grid  capability  are 
used  in  conjunction  with  the  robust  SIMPLE  algorithm  and  the  enthalpy  formulation  to 
faithfully  simulate  this  complex  process. 


CHAPTER  2 

FORMULATION  FOR  CONVECTIVE  HEAT  TRANSFER 
AND  FREE  SURFACE  FORMATION 

2. 1  Background 

In  this  chapter,  a  formulation  is  developed  for  convection  heat  transfer  and  free  surface 
formation  in  float  zones.  Of  particular  interest  is  natural  convection  generated  by  density 
differences  in  the  fluid,  and  thermocapillary  convection  caused  by  shear  stress  arising  from 
thermal  and  composition  gradients  along  a  free  surface  separating  a  liquid  and  a  gas.  The 
interest  in  such  flows  arises  from  its  prevalence  in  nature  as  well  as  in  engineering 
applications  and  devices.  For  instance,  natural  convection  is  one  of  the  important  and 
commonly  employed  method  to  cool  electronic  devices.  Thermocapillary  convection 
becomes  important  in  many  crystal  growth  and  solidification  processes,  notably,  the  float 
zone  process  for  single  crystal  growth  from  the  melt.  The  chapter  starts  with  a  description 
of  the  physical  scales  that  characterize  natural  and  thermocapillary  convection,  illustrating 
the  nondimensional  parameters  that  determine  the  strength  and  relative  importance  of  these 
physical  processes.  Then  the  governing  equations  are  transformed  into  a  body-fitted 
coordinate  system  which  is  convenient  for  handling  curved  and  complex  geometry.  In  this 
work,  the  governing  equations  are  solved  by  means  of  the  popular  and  robust  SIMPLE 
algorithm  for  incompressible  flows.  For  convection  dominated  flows,  the  numerical 
treatment  of  the  convection  terms  in  the  Navier-Stokes  equations  has  a  substantial  impact 
on  the  robustness  and  accuracy  as  will  be  demonstrated  in  the  next  chapter.  This  chapter 
concludes  with  a  formulation  for  the  generation  of  meniscus  shape  for  float  zone 
configurations.  The  issues  of  existence,  multiplicity  and  sensitivity  of  free  surface  shapes 
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are  discussed  in  the  context  of  solutions  to  the  Young-Laplace  equations  and  a  free  energy 
concept  originated  by  Shyy  et.  al.  (1993). 

2.2  Sources  of  Convection 

2.2. 1    Natural  Convection 

The  strength  of  natural  convection  depends  on  the  density  differences  caused  by  thermal 
and  compositional  inhomogeneities  within  the  fluid.  In  two-dimensions,  the  source  of 
buoyancy  driven  convection  is  the  component  of  the  local  density  gradient  normal  to  the 
gravity  vector.  A  temperature  scale,  AT,  can  be  defined  which,  along  with  a  characteristic 
length  scale,  d,  characterizes  the  horizontal  thermal  gradient  and  consequently  the  horizontal 
density  gradient.  In  addition,  we  can  define  a  characteristic  solute  scale,  Ac|).  The  Boussinesq 
approximation  is  now  invoked,  which  by  linearization  neglects  the  density  variation  in  all 
the  terms  of  the  governing  equations  except  the  source  term.  The  buoyancy  source  term  then 
becomes 

"  (Q  "  Qref)  g  =  0  g  [  P  T  (T  "  Tref)  +       (<>  -  cj)ref)]  (2.1) 

where  |3T  and  are  the  coefficients  of  thermal  expansion  and  solutal  expansion, 
respectively,  of  the  melt.  The  resulting  nondimensional  parameters  are 

Thermal  Grashof  number: 

Solutal  Grashof  number: 
Prandtl  number: 
Schmidt  number: 

where  k  is  the  thermal  conductivity,  q  is  the  density,  Cp  is  the  heat  capacity,  \i  is  the  dynamic 
viscosity,  v  is  kinematic  viscosity,  and  D  is  the  mass  diffusivity.  The  thermal  and  the  solutal 
Rayleigh  numbers  can  be  defined  as  RaT  =  GrT  •  Pr  and  Ra^  =  Gr^  •  Sc,  respectively. 


GrT  = 


g  PT  AT  d3 


Gr«,  = 


g  Pj>  Acj>  d3 


Pr  = 


k/q 


(2.2a) 


(2.2b) 
(2.2c) 
(2.2d) 
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2.2.2    Marangoni  Convection 

For  problems  involving  free  or  deformable  interfaces  between  two  liquids  or  between 
a  liquid  and  a  gas/vapor,  temperature  and  composition  inhomogeneities  can  result  in  a 
surface  tension  gradient  along  the  interface.  This  exerts  a  shear  stress  on  the  fluid  giving  rise 
to  thermocapillary  convection  also  known  as  Marangoni  convection.  The  float  zone  is  an 
example  of  a  configuration  where  Marangoni  convection  plays  an  important  role  in 
governing  the  solidification  characteristics  of  the  melt.  General  discussion  of  Marangoni 
convection  can  be  found  in  Koschmieder  (1993),  Shyy  (1994).  We  now  briefly  examine  the 
case  of  thermally  generated  Marangoni  convection. 

For  liquids  bounded  by  a  free  or  deformable  boundary  where  a  temperature  gradient 
exists  along  the  interface,  a  surface  shear  stress  is  generated  by  the  surface  tension  gradient. 
Thus,  the  force  balance  can  be  written  as 

\i  ■  V  (V  •  T)  •  n  =  |£  VT  •  T  (2.3) 

where  the  surface  tension  gradient  is  proportional  to  the  temperature  gradient  and  gives  rise 
to  the  normal  derivative  of  the  tangential  velocity  at  the  free  surface,  y  is  the  surface  tension, 

V  is  the  velocity  vector,  n  is  the  unit  normal  vector  to  the  free  surface,  and  T  is  the  unit 

tangent  vector  to  the  free  surface.  The  convection  strength  is  described,  in  nondimensional 
terms,  by  the  Marangoni  number, 


Ma  = 


dy 


dT 


ATd 

  (2.4) 


where  a  is  the  thermal  diffusivity.  The  Marangoni  number  represents  the  balance  between 
the  shear  stress  and  the  surface  tension  gradient  after  the  velocity  and  space  variables  have 
been  nondimensionalized  by  the  velocity  scale  based  on  heat  conduction,  as  will  be 
described  in  the  next  section. 
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2.3  Nondimensionalization  and  Scaling  Procedure 
The  governing  equations  can  be  nondimensionalized  with  reference  length,  velocity  and 
thermal  scales.  The  choice  of  nondimensional  parameters  used  to  normalize  the  governing 
equations  is  not  unique  as  will  be  apparent  in  the  following. 
The  dimensionless  variables  can  be  defined  as: 

t  =  t/tr  ,  X  =  x/d  ,  Y  =  y/d 

u  =  u/Ur  ,  v  =  v/Ur  (2.5) 
P  =  p/(QUr2)  and  0  =  (T  -  Tc)/(Th  -  Tc) 

In  equation  (2.5),  d  is  the  reference  length  scale,  Ur  is  the  reference  velocity  scale,  Tf,  and 
Tc  are  the  maximum  and  minimum  temperatures  present  in  system,  and  the  temperature  is 
normalized  by  the  reference  temperature  scale,  AT  =  Th  —  Tc,  such  that  the 
nondimensional  temperature  is  between  zero  and  unity.  The  reference  time  scale  is  defined 
astr  =  d/Ur.  The  reference  velocity  scale  is  determined  by  the  dominant  physical  processes 
occurring  in  the  system  as  described  next. 

2.3. 1  Heat  Conduction  Scales 

The  reference  velocity  scale  for  heat  conduction  is 

Ur  =  §  (2.6) 

where  a  =  k/ (qCp)  is  the  thermal  diffusivity  of  the  fluid.  Then,  with  the  reference  scales 
specified  earlier,  the  non-dimensional  equations  become: 

momentum:  ^  +  V  ■  V(V)  =  -  VP  +  PrV2V  -  RaPr0g  (2.7a) 

energy:  ^  +  V  •  V0  =  V20  (2.7b) 

where  V  is  the  nondimensional  velocity  vector,  Pr  is  the  Prandtl  number,  g  is  the  unit  vector 
in  the  direction  of  gravity,  and  Ra  is  the  thermal  Rayleigh  number. 

2.3.2  Natural  Convection  Scales 

The  reference  velocity  scale  for  natural  convection  is 
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Ur  =  yg  (3  d  AT  (2.8) 
and  with  the  reference  scales  specified  earlier,  the  nondimensional  equations  become 


momentum: 


dV 

dx 


+  v  •  V(V)  =  -  VP  +  -^=V2V  -  0g  (2.9a) 

7GrT 


energy:  +  y  •  V0  =  j=  V20  (2.9b) 

*  PrTGr^ 

where  Pr  is  the  Prandtl  number  and  Grj  is  the  thermal  Grashof  number. 

2.3.3  The  Marangoni  Number 

If  equation  (2.3)  is  nondimensionalized  with  respect  to  the  heat  conduction  scales  (2.6), 
the  Marangoni  number  defined  by  equation  (2.4)  occurs  naturally  in  the  dimensionless  form 

V  (V  •  T)  •  n  =  Ma  V0  ■  T  (2.10) 

Considering  that  the  dimensionless  temperature  gradient  is  of  order  one,  the  above  equation 
indicates  that  the  strength  of  the  shear  stress  in  the  free  surface  region  is  proportional  to  the 
Marangoni  number;  this  shear  stress,  in  turn,  induces  convection  within  the  body  of  the  fluid. 
Thus,  the  Marangoni  number  describes  the  strength  of  convection  due  to  the  surface  tension 
gradient. 

2.3.4  Marangoni  Convection  Scales 

For  a  flow  where  viscous  effects  dominate,  the  reference  velocity  scale  for 
thermocapillary  convection  is 

(dv/dT)  AT 

U~^y   (2.11) 

whereas  for  a  highly  convective  flow, 
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2.4  Governing  Equations  and  Solution  Procedure 
The  governing  equations  in  Cartesian  coordinates  for  two-dimensional,  incompressible 
flow  can  be  written  in  dimensional  form  as: 


continuity: 


x-momentum: 


d(Qu)     d(Quu)  d(puv) 


at 


ax 


ay 


dg  a(Qu)  a(gv) 
at     ax  ay 


o 


ax 


1L  iu\+ iu 
axl^  ax/    ay  I  ^  ay 


+  Su(x,y) 


(2.13a) 


(2.13b) 


v-momentum: 


a(Qv)    a(Quv)  a(Qw) 


at 


ax 


ay 


3y 


-  (Q  -  Qref)  §  +  Sv(x,y) 


ay, 


(2.13c) 


energy: 


a(gCpT)    a(QuCpT)  a(QvCpT) 


at 


ax 


9y 


axlK  axj  +  ay\K  ayj 


+  ST(x,y) 


(2.13d) 


In  the  above,  the  dimensional  form  has  been  chosen  because  it  is  convenient  to  represent 
variable  physical  properties  such  as  q,  the  density,  u,  the  viscosity,  k,  the  conductivity,  and 
Cp,  the  specific  heat  in  the  conservative  form.  Here,  u  is  the  velocity  component  in  the 
x-direction,  v  is  the  velocity  component  in  the  y-direction,  p  is  the  pressure,  T  is  the 

A 

temperature,  and  g  is  acceleration  due  to  gravity,  The  gravity  vector  g  =  -  g  j  is  assumed 
to  be  in  the  -y  direction.  For  inclined  configurations,  the  appropriate  components  should 
appear  in  both  the  momentum  equations.  Su,  Sv  and  St  are  appropriate  source  terms  in  the 
governing  equations.  The  form  of  the  source  terms  is  problem  dependent  and  will  be 
described  later  in  the  appropriate  context. 
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It  frequently  is  convenient  to  write  and  solve  the  governing  equations  in  generalized, 
nonorthogonal,  curvilinear  coordinates,  especially  to  fit  complex  geometries  and  implement 
the  boundary  conditions  on  surfaces  or  walls  with  arbitrary  orientations.  Furthermore,  in 
many  applications  it  is  desirable  to  preferentially  improve  grid  distribution  in  certain  regions 
of  the  flow.  These  regions  may  be  determined  manually  or  may  evolve  dynamically  during 
the  solution  process. 
The  transformation  may  be  expressed  as 


I  =  l(x,y) 
r\  =  T](x,y) 

The  continuity  equation  in  nonorthogonal  coordinates  then  becomes, 

^(Je)  +  ^H)  +  ^v)  =  0 


(2.14a) 


(2.14b) 


where 

J  =  -  X^ 

is  the  Jacobian  of  the  transformation  and 

U  =  uyr,  -  vxn 
V  =  vx^  —  uy^ 


(2.14c) 


(2.14d) 


are  the  contravariant  velocity  components.  The  momentum  equations  then  become 

y(qn^  -  q2UT)) 


d(Jou)  +  d(oUu)  +  a(QVu) 


at 


91 


dp  dp 


3? 


y(q2u?  -  qsuTi) 


+  SU(1,T|)  •  J 


(2.14e) 


a(Jov)  +  a(oUv)  +  a(oVv) 


at 


dp  dp 

xsa^  ~  Xt1af 


+  4 
a^ 


+  A 

dr\ 


^(q,v?  -  q^) 
j(q2v§  -  q3vn)    -  (C  -  Qref)  E  '  J  +  Sv  (1,  Tl)  •  J 


(2.14f) 


along  x  and  y  directions  respectively,  where  x^,  y^  and  so  on  are  the  metrics  of  the 
transformation  and  qi,  q2  and  q3  are  defined  as 
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qi  =  A  +  y\ 

q2  =  x^Xt,  +  y^yr,  (2.14g) 
q3  =  x|  +  y| 

The  energy  and  other  scalar  transport  equations  can  be  written  down  in  a  similar  manner  and 
will  not  be  repeated  here.  The  advantage  of  this  form  is  that  complex  geometries  can  be  fitted 
by  body-conforming  grids,  whereby  boundary  conditions  can  be  applied  easily  and  the  grid 
can  be  distributed  preferentially  in  desired  regions  so  as  to  obtain  better  resolution. 

A  pressure-based  numerical  procedure,  SIMPLE,  has  been  chosen  for  computing 
laminar  incompressible  flows  of  interest  in  this  work.  The  details  of  the  basic  pressure 
correction  algorithm  are  given  in  Patankar  (1980),  with  the  extension  of  the  procedure  to 
general  curvilinear  coordinates  given  in  Shyy  (1994).  The  present  implementation  of  the 
algorithm  follows  the  work  of  Braaten  and  Shyy  (1986)  and  Shyy  et  al.  (1985).  A  staggered 
grid  arrangement  is  adopted  for  the  discretization  of  Eqs.  (2.14b-f)  with  the  layout  of  the 
flow  variables  illustrated  in  Fig.  2.1.  The  finite  volume  form  of  the  conservation  laws  may 
be  obtained  by  integrating  Eqs.  (2. 14b-f)  over  the  sample  control  volume  shown  in  Fig.  2. 1 . 
The  diffusion  terms  are  always  approximated  using  second  order  central  differences  in  this 
study  and  the  time  derivative  by  a  first  order  backward  Euler  scheme.  A  variety  of  convection 
schemes  have  been  used  and  a  detailed  discussion  is  provided  later  in  the  next  section.  In  the 
SIMPLE  algorithm,  a  segregated  approach  is  adopted  wherein  the  momentum,  energy  and 
other  scalar  equations  are  solved  in  a  sequential  manner.  A  pressure  correction  approach 
(Patankar  1980,  Shyy  1994)  is  used  to  enforce  mass  continuity  in  each  discrete  control 
volume  and  the  velocity  components  are  corrected  by  the  consistent  and  efficient  D'yakonov 
procedure  discussed  in  detail  in  Shyy  (1994). 

2.5   Formulation  and  Methodology  for  Free  Surface  Computation 

Now,  it  is  appropriate  to  take  up  the  study  of  the  generation  of  meniscus  shapes  for  float 
zone  configurations.  The  formation  of  the  meniscus  during  the  melting  process  is  a  major 
concern  in  the  successful  operation  of  the  float  zone  technique  and  this  aspect  will  be  the 
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Figure  2.1.  Layout  of  the  dependent  variables  on  the  staggered  grid.  The 
direction  of  the  contravariant  velocities  are  different  from  the  cartesian 
components.  The  contravariant  components  are  perpendicular  to  the  respective 
grid  lines  on  which  they  are  defined. 

considered  in  detail  in  the  next  chapter  by  means  of  several  examples.  Specifically,  the  issues 
of  existence,  multiplicity  and  sensitivity  of  solutions  to  the  Young-Laplace  equation 
representing  the  liquid  meniscus  formed  under  equilibrium  conditions,  subject  to  different 
combinations  of  Bond  number,  aspect  ratio  and  pressurization,  will  be  addressed  in  this 
chapter  and  investigated  by  mean  of  several  examples  described  in  the  next  chapter.  A  free 
energy  concept  originated  by  Shyy  et.  al.  (1993)  is  described  next  and  will  be  applied  to 
address  the  physical  realizability  of  a  given  meniscus  shape.  The  approach  is  static  in  nature, 
which  is  a  satisfactory  approximation  considering  the  slow  speed  of  the  movement  of  the 
melt.  Fig.  1 . 1  is  a  schematic  illustration  of  the  float  zone  process  and  Fig.  2.2  is  a  schematic 
of  the  mathematical  model,  which  allows,  as  it  should,  the  trijunction  point  to  move  to  satisfy 
a  given  value  of  the  static  contact  angle  or  other  constraints  such  as  a  fixed  volume  of  the 
melt. 

Meniscus  profiles  obey  the  Young-Laplace  equation  in  axisymmetric  form  (Finn  1986, 
Myshkis  et.  al.  1987,  Shyy  1994),  which  relates  the  curvature  of  the  free  surface  to  the 
pressure  difference  across  it.  The  governing  equation  describes  the  generation  of  surface 
curvature  to  balance  the  sum  of  the  hydrostatic  pressure  and  pressurization.  Specialized  to 
an  axisymmetric  geometry,  it  can  be  written  as  follows: 
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Figure  2.2.  Schematic  of  the  mathematical  model.  Major  simplifications 
involve  isothermal  melt  zone  and  a  flat  solid-liquid  interface  as  opposed  to  the 
curved  one  shown  in  Fig.  1.1.  Axisymmetry  has  been  assumed. 


where  Aq  is  the  density  difference  between  the  ambient  and  the  melt  and  AP  is  the 
pressurization.  y  is  the  surface  tension  coefficient  and  the  meniscus  is  described  by  f  =  r(y), 
where  r  is  the  radius  and  y  is  the  vertical  coordinate.  A  variety  of  boundary  conditions  have 
been  treated  in  literature,  such  as  fixed  contact  angles  at  both  top  and  bottom  (Finn  1986, 
Concus  1990,  Myskis  et.  al.  1987),  fixed  radius  at  the  bottom  and  fixed  contact  angle  at  the 
top  (Shyy  et.  al.  1993)  and  fixed  radius  at  both  top  and  bottom  (Coriell  et.  al.  1977,  Li  et. 
al.  1993,  Kobayashi  1984,  Kobayashi  1988,  Zhang  and  Alexander  1992,  Murthy  1987). 
These  conditions  have  been  devised  for  configurations  under  different  physical  constraints. 
In  the  context  of  floating  zones,  the  boundary  conditions  to  be  imposed  may  vary  depending 
on  the  operating  conditions.  Therefore  this  study  has  been  conducted  with  a  variety  of 
different  constraints  imposed  on  the  solutions.  The  following  boundary  conditions  may  be 
imposed: 


Aq  g  y  -  AP  =  y 


r  i 


(2.15a) 


(1  +  f'2)3/2      f  (1  +  f'2)l/2 


f(y=0)  =rb 


(2.15b) 


either 


f  (y=hc)  =  rc 


(2.15c) 
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or  f  (y=hc)  =  tan(  ji/2  -  (J)c)  (2.15d) 

Figure  2.2  shows  a  schematic  of  the  geometry.  For  an  isothermal  liquid  bridge,  Eq.  (2. 15a) 
can  be  nondimensionalized  with  a  length  scale  of  rb  and:a  pressure  scale  of  y/r^,  to  yield, 


Bo  y  -  AP  = 


1 


'2  '2 
.(1  +  f  )V2      f  (1  +  f  )l/2_ 


(2.16) 


where  the  italics  indicate  non-dimensionalized  quantities.  Bo  is  the  Bond  number  defined 

as: 

Aq  g  rj 

Bo  =  — y-2  (2.17) 

In  Shyy  et.  al.  (1993),  the  nonlinear  differential  equation,  Eq.  (2.16),  has  been  solved 
subject  to  a  fixed  radius  at  the  lower  boundary  and  a  fixed  contact  angle  at  the  top  boundary. 
Such  a  two  point  boundary  value  problem  cannot  yield  an  unique  solution.  Thus  there  arises 
the  problem  of  selecting,  out  of  the  multiple  solutions  that  are  mathematically  permissible, 
the  one  that  corresponds  to  the  physical  equilibrium  condition  that  exists  in  reality. 
Therefore  we  now  invoke  the  thermodynamic  condition  that  at  equilibrium  the  free  energy 
is  a  minimum  and  select  the  solution  that  minimizes  the  free  energy. 

For  simplicity,  consider  only  the  isothermal  condition  with  no  phase  change  and 
examine  the  Helmholtz  free  energy,  F  =  U  -  TS  which  contains  three  contributions: 
(i)  The  potential  energy  from  the  effective  head 


Uj  =  |jrfz  (Aggy  -  AP)  dy  (2.18a) 

(ii)  The  surface  energy  of  the  meniscus  forming  the  gas-liquid  interface 

U2  =  J  2ji  Y  f  (l  +  f'2)1/2  dy  (2.18b) 

(iii)  The  surface  energy  needed  to  wet  the  solid-liquid  wetted  area 

u3  =  (Yis  -  Ysg)  n  r2.  (2.18c) 
and  from  Young's  contact  condition  this  becomes, 
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U3  =  -  y  cos(J)0  n  xl  (2.18d) 
where  (j)0  is  the  static  contact  angle.  The  equilibrium  shape  is  the  one  that  minimizes  the 
total  free  energy: 

he 

E  =  |  [  n  f 2  (Afjgy  -  AP)  +  2ft  Y  f  (1  +  f'2)1/2]  dy  -  Jt  y  I"2  cos(())0)  (2.18e) 
o 

The  Young-Laplace  equation  may  be  integrated  starting  with  a  fixed  radius  and  base  angles, 
—  Jt/2  <  cpb  <  jt/2,  to  yield  all  possible  solutions  for  a  given  aspect  ratio  of  the  domain. 
Subsequently,  solutions  satisfying  a  given  contact  angle  at  the  top,  (J)c,  are  selected  and  their 
free  energies  calculated  according  to  Eq.  (2. 18e).  If  any  of  these  solutions  locally  minimize 
the  free  energy,  then  these  profile  shapes  are  statically  stable. 

As  described  in  Shyy  et.  al.  (1993)  and  Shyy  (1994),  the  value  of  (j)0  is  specified,  and 
E  is  calculated  by  fixing  the  value  of  <}>0  in  Eq.  18e)  and  scanning  through  the  whole  range 
of  meniscus  profiles  obtained  by  fixing  the  lower  trijunction  location  (with  varying  angle) 
and  the  height  of  the  upper  trijunction  point  (with  varying  locations  and  angles),  distinct 
extrema  are  obtained  at  specified  values  of  §0  .  Among  the  multiple  solutions,  the  one 
corresponding  to  a  point  of  minimum  on  the  curve  is  the  physically  realizable  stable 
meniscus  profile.  When  either  a  maximum  or  a  non-extremum  point  arises,  the  solution  is 
unstable. 

The  volume  of  the  melt  zone  in  non-dimensional  form  is  given  by, 

[  ^  A 

V  =       jt  f  dy  (2.19a) 


o 

where 

=  hc/rb  (2.19b) 
is  the  aspect  ratio  of  the  zone.  In  float  zones,  the  meniscus  profiles  may  be  required  to  enclose 
a  fixed  volume  of  the  melt  thus  satisfying  Eq.  (2. 19)  (Coriell  et.  al.  1977,  Finn  1986,  Myshkis 
1987,  Zhang  and  Alexander  1992,  Li  et.  al.  1993).  Furthermore,  non-isothermal  liquid 
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bridges  need  to  be  considered  considered  for  the  float  zone,  see  for  instance  Zhang  and 
Alexander  (1992).  However,  the  present  study  described  in  the  next  chapter,  will  be 
restricted  to  the  isothermal  case  to  simplify  the  number  of  parameters  involved.  The  aim  is 
to  determine  the  sensitivity  of  the  profile  shapes  to  the  imposed  boundary  conditions.  As 
will  be  demonstrated,  the  basic  considerations  from  free  energy  and  uniqueness  can  shed 
useful  light  on  the  issues  relevant  to  the  float  zone.  It  is  also  noted  that  additional  forces  may 
be  generated  on  the  free  surface  due  to  electromagnetic  effects  if  the  zone  is  heated  by  an 
induction  coil  (Lie  et.  al.  1989,  Riahi  and  Walker  1989). 

Dynamic  effects  have  been  considered  as  secondary  effects  for  the  purposes  of  this 
analysis;  however,  they  are  incorporated  later  in  this  work  with  fully  coupled  moving  grid 
calculations  with  fluid  convection,  phase  change  and  free  surface  dynamics.  Dynamic 
effects  have  been  also  considered  by  several  investigators  (Langbein  1992,  Tsamopoulos  et. 
al.  1988  and  Tsamopoulos  et.  al.  1992).  Quasiequilibrium  studies  involving  contact  angle 
hysterisis  have  been  treated  in  Shyy  et.  al.  (1993)  and  Shyy  (1994).  It  is  considered  that  the 
phenomenon  of  contact  angle  hysterisis  is  largely  due  to  surface  roughness  (Huh  and  Mason 
1977)  and  surface  inhomogeneity  (Shyy  1994).  In  the  following,  it  is  considered  that  the 
melt  zone  moves  slowly  and  a  static  consideration  suffices. 

2.6  Results  and  Discussion 

Meniscus  profiles  were  obtained  by  solving  the  Young-Laplace  equation  in  cylindrical 
coordinates  subject  to  the  boundary  conditions  of  fixed  radius  at  the  bottom  and  either  (I) 
fixed  contact  angle  at  the  top,  or  (II)  fixed  volume  of  the  liquid  bridge  (half  domain  assuming 
axial  symmetry). 
2.6. 1    Low  Bond  Number 

Initial  calculations  were  conducted  following  the  model  proposed  by  Shyy  et.  al.  (1993) 
and  compared  with  experimental  data  provided  by  Dreeben  et.  al  (1980).  The  parameters 
were  as  follows:  rb  =  2.42  X  1(H  m,  Bo  =  2.825  x  10~3,  AR  =  2.24,  contact  angle,  <{>c  =115 
deg  and  pressurization,  AP  =  4.35%  of  the  ambient  pressure  given  by  Pambient  =  1.0  x  105 
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Nnr2.  In  this  chapter,  AP  is  always  specified  as  a  percentage  of  Pambient  =  1-0  x  105  Nm-2. 
Here  rb  is  the  radius  of  the  feed  rod  and  it  is  assumed  that  the  meniscus  profile  is  pinned  to 
the  edge  of  the  feed  rod.  It  may  be  noted  from  the  photographs  of  floating  zones  provided 
in  Dreeben  et.  al.  (1980)  that  the  contact  angle  is  not  constant  over  the  contact  line  nor  is  the 
liquid  bridge  axisymmetric.  The  solid-liquid  interface  is  curved  and  not  flat  as  is  assumed 
in  the  simplified  model  used  in  the  present  calculations.  It  may  also  be  noted  that  the 
variation  of  surface  tension  with  temperature  is  also  not  considered  in  this  model.  The 
material  under  consideration  is  Sapphire  (AI2O3)  whose  properties  are  listed  in  Table  2.1. 

Figure  2.3  shows  the  meniscus  profiles  obtained  using  the  theoretical  model  as  well  as 
the  experimental  data.  In  Fig.  2.3,  (i)  shows  the  actual  profiles  obtained  by  integrating  the 
Young-Laplace  equation,  that  satisfy  the  imposed  contact  angle  condition  at  the  top.  (ii) 
shows  the  rc  versus  (J)c,  (iii)  shows  (})c  versus  (j)b  and  (iv)  is  the  free  energy  curve,  E  versus 
(j)c.  In  (i),  the  experimental  data  from  Dreeben  et.  al.  (1980)  is  shown  alongside  for 
comparison  purposes,  (ii)  and  (iii)  show  the  range  of  solutions  which  may  be  obtained  from 
the  Young-Laplace  equation.  It  is  obvious  that  there  are  two  solutions  satisfying  the  given 
contact  angle  condition.  The  profile  that  minimizes  the  free  energy  is  shown  with  a  solid 
line  and  it  matches  the  experimentally  observed  profile  closely. 

A  sensitivity  study  has  been  conducted  to  study  the  formation  and  physical  realizability 
of  the  meniscus  profiles  subject  to  various  boundary  conditions  and  parameters.  First, 
holding  the  Bond  number  constant,  the  aspect  ratio  was  varied  over  the  range 
0.5  <  /R  <  4.0. 

Figure  2.4  shows  the  profile  shapes  at  each  aspect  ratio.  Figure  2.5  is  a  summary  of  the 
results  obtained  by  imposing  a  desired  contact  angle,  (j)c  =  115  deg  for  the  various  aspect 
ratios  considered.  Multiple  solutions  may  exist  especially  at  lower  aspect  ratios  The 
available  solutions  are  marked  with  an  '  x '  and  the  solutions  having  the  lower  free  energy 
are  marked  with '  <g) ' .  A  notable  characteristic  is  the  rapid  variation  in  the  meniscus  profiles 
for  aspect  ratios  in  the  range,  2.25  <  /R  <  2.5.  The  value  of  rc  which  shows  a  decreasing 
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Bo  =  2.825  X  10" 3  ,  AP  =  4.35%  ,  Pambient  =  1.0  X  105Nm-2  ,        =  2.24 


  Computation 

  Computation 

o-  -  -  ©  Experiment 


Figure  2.3.  Comparison  of  calculated  meniscus  shapes  with  the  experimental  data  of 
Dreeben  et.  al.  (1980).  (i)  Calculated  and  experimentally  obtained  profiles.  It  may  be 
observed  that  two  solutions  are  obtained  with  the  given  contact  angle  of  1 15  deg. 
The  profile  with  the  lower  free  energy  is  closer  to  the  experimentally  observed  profile, 
(ii)  rc/rb  versus  (pc  (iii)  cpc  vs  (pb  (iv)  Free  energy,  E  vs  (pc 


trend  with  increasing  aspect  ratio  up  to  M  =  2.25  shows  a  sudden  increase  as  the  aspect  ratio 
increases  to  2.5.  The  decreasing  trend  resumes  thereafter.  Only  the  profiles  having  the  lower 
free  energy  are  shown.  Profiles  with  aspect  ratio  of  2.5  and  higher  show  a  'neck'  or  turning 
point. 
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<|)c  =  1 15  deg,  Bo  =  2.825  x  10 


3.75  ; 

3.5  : 
3.25  1 


1.75  : 

i  -5  : 

1.25  ■ 

1  ■ 

0.75  • 

0.5  ': 

0.25  : 


Figure  2.4.  Profiles  satisfying  contact  angle  of  115  deg.  for  Bo  =  2.825  x  10~3  for 
0.5  <  /R  <  4.0  Only  profiles  with  the  lower  free  energy  are  shown. 


Table  2.1  Bond  numbers  for  materials  of  interest,    =  1.25  x  10~2  m. 
Data  obtained  from  Keene  (1993)  and  Li  et.  al.  (1993). 


material 

density 

surface  tension 

Bond  Number 

A1203 

3400  kg  nr3 

0.69  N  nr1 

7.55 

GaAs 

5720  kg  m-3 

0.455  N  nr1 
(  w.r.t  B203 ) 

13.72 
(with  encapsulant) 

NiAl 

5000  kg  nr3 

1.35  Nnr1 

5.68 
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rc/rb  vs.  AR 


Constant  <|>  =  1 15° 


Bo  =  2.825  x  10" 


V  vs.  AR 


Figure  2.5.      Solutions  for  an  imposed  contact  angle  of  1 15  deg  in  the  range 
0.5  <  /R  <  4.0 

(i)  rc  vs.  the  aspect  ratio.  Note  the  change  in  the  trend  of  the  solutions  when  the  aspect 
ratio  changes  from  2.25  to  2.5.  (ii)  (J>b  that  give  rise  to  the  imposed  contact  angle, 
(iii)  volume  vs.  the  aspect  ratio.  In  all  cases,  if  multiple  solutions  are  obtained,  then 
the  solutions  with  the  lower  free  energy  are  marked  with  a  '  <8> .' 


With  a  set  pull  rate  and  heater  power,  at  steady  state,  the  float  zone  occupies  a  prescribed 
volume.  Hence  it  is  useful  to  examine  the  solution  characteristics  subject  to  the  constraint 
of  a  fixed  volume.  However,  in  reality,  the  change  in  heat  transfer  characteristics  due  to  the 
coupling  between  the  movement  of  the  float  zone  and  heat  transfer  due  to  conduction  as  well 
as  buoyancy  driven  flow,  may  result  in  changes  of  the  prescribed  melt  volume  with  time. 
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Furthermore,  the  issue  of  viewing  the  static  contact  angle  as  a  material  property  needs  to  be 
addressed  also. 

Figure  2.6  displays  free  energy  curves  at  selected  aspect  ratios.  It  may  be  seen  that  all 
the  selected  curves,  except  the  case  of  /R  =  1.75,  show  a  local  minimum  at  (J)0  =  115  deg. 
The  profile  shown  at  /R  =  1 .75  cannot  achieve  static  stability  for  the  specified  contact  angle. 


Figure  2.6.  Energy  versus  (j)c  curves  for  selected  aspect  ratios  and  4>0  =  115  deg. 
It  may  be  noted  that  at  an  aspect  ratio  of  1.75,  the  free  energy  curve  shows  a  local 
maximum  and  hence  cannot  be  statically  stable.  All  other  curves  display  a  local 
minimum  indicating  static  stability. 
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It  is  interesting  to  observe  that  such  a  "window"  of  exists  within  which,  although  the 
mensci  exist  mathematically,  they  are  not  physically  realizable.  Above  and  below  this 
"window"  statically  stable  menisci  exist.  In  particular,  it  may  be  noted  that  profiles  at  /R  = 
2.5  and  2.75  exhibit  'necking.'  The  stability  implication  of  such  a  shape  needs  to  be 
addressed  in  a  dynamic  framework.  Here,  the  necking  can  satisfy  the  static  stability 
requirement. 

Figures  2.7  and  2.8  are  the  corresponding  solutions  that  satisfy  a  fixed  volume  criterion 
of  V  =  5.908,  chosen  so  as  to  correspond  with  the  theoretical  solution  shown  in  Figure  2.3. 
In  this  case  the  solution  characteristics  vary  smoothly  with  the  aspect  ratio.  However,  the 
rc  shows  a  qualitatively  different  trend  as  the  aspect  ratio  changes  from  2.25  to  2.5.  The 
corresponding  meniscus  profiles  shown  in  Figure  2.8  start  to  exhibit  a  'neck'  in  the  profile 
for  aspect  ratios  of  2.5  and  above,  although  the  'necks'  are  shallower  than  those  seen  in  the 
fixed  contact  angle,  (J)c  =  115  deg  case  shown  in  Figure  2.4.  The  contact  angle,  4>c, 
correspondingly  exhibits  a  sharp  decrease  as  the  aspect  ratio  increases  from  2.25  to  2.5.  This 
is  consistent  with  the  appearance  of  turning  points  or  'neck'  in  the  profile  shapes  as  shown 
in  Figure  2.8. 

The  next  case  shown  in  Figure  2.9  examines  the  sensitivity  of  the  solutions  to  changes 
in  the  Bond  number.  The  Bond  number  was  doubled  to  Bo  =  5.650  x  10-3  retaining  the 
pressurization  of  4.35%.  Figure  2.9  shows  the  range  of  the  available  solutions.  It  is  readily 
noted  that  there  is  a  substantial  reduction  in  the  range  of  the  available  solutions.  The 
Young-Laplace  equation  cannot  be  integrated  for  aspect  ratios  greater  than  0.75.  The  (f)c 
versus  (j)b  curves  show  that  multiple  solutions  do  not  exist  for  the  aspect  ratios  considered. 
It  may  be  concluded  that  for  the  parameter  ranges  and  geometrical  configurations  considered 
that  the  meniscus  profiles  are  most  sensitive  to  variations  in  Bond  number  rather  than 
variations  in  the  aspect  ratio  or  the  different  boundary  conditions  used  in  the  calculations. 
This  set  of  computations  can  help  explain  the  non-existence  of  floating  zones  as  the  surface 
tension  coefficient  is  not  sufficient  to  counteract  the  hydrostatic  pressure  of  the  melt. 
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V  =  5.908,  Bo  =  2.825  x  10 


Figure  2.7  Profiles  satisfying  constant  volume,  V  =  5.908,  for  Bo  =  2.825  x  10-3 
for  0.5  <  /R  <  4.0  Only  profiles  with  the  lower  free  energy  are  shown. 


2.6.2    High  Bond  Number 

In  practice,  the  geometrical  configurations  employed  result  in  a  high  value  of  the  Bond 
number.  In  particular,  we  are  interested  in  the  following  configuration.  Given  a  feed  rod 
diameter  of  25  mm,  Table  2. 1  shows  the  surface  tension  coefficients  and  the  corresponding 
Bond  numbers  for  some  materials  of  current  interest. 
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Figure  2.8      Solutions  for  an  fixed  volume  of  5.9  in  the  range  0.5  <  /R  <  4.0 

(i)  rc  versus  the  aspect  ratio 

(ii)  (j)c  versus  the  aspect  ratio 

Note  the  change  in  the  trend  of  the  solutions  when  the  aspect  ratio  changes  from 
2.25  to  2.5.  (iii)  (J>b  that  give  rise  to  solutions  satisfying  V  =  5.9.  In  all  cases, 
if  multiple  solutions  are  obtained,  then  the  solutions  with  the  lower  free  energy 
are  marked  with  a '  <S> '  . 


2.6.3    Determination  of  Pressurization  Values 

The  meniscus  shape  is  substantially  influenced  by  the  value  of  AP,  pressurization,  in  Eq. 
(2.15).  In  literature,  the  value  of  the  pressurization  is  usually  determined  by  satisfying  a 
constant  volume  criterion  in  addition  to  the  boundary  conditions  imposed.  In  the  present 
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Figure  2.9.  Solutions  for  Bo  =  5.65  x  10-3.  Solutions  are  obtained  for  aspect 
ratio  of  0.5  and  0.75.  The  decrease  in  the  number  of  available  solutions  may  be 
noted  as  opposed  to  the  case  of  Bo  =  2.825  x  10-3. 

(i)  (j)c  vs.  (j>b 

(ii)  rc  vs.  (})c  and 

(iii)  volume  vs.  (j)c. 


calculations,  the  objective  was  to,  instead,  find  a  range  of  pressurization  values  for  which 
aspect  ratios  close  to  unity  can  be  achieved.  The  result  of  such  a  search  is  described  in  Fig. 
2.10.  A  Bond  number  of  8.8  was  selected  in  order  to  effect  a  comparison  with  the  results 
of  Quon  et.  al.  (1993).  The  pressurization  ranges  obtained  from  Fig.  2.10  are  used  to 
generate  the  subsequent  solutions. 
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Bo  =  8.8 

1.4 
1.2 
1 

10.8 
<  0.6 
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0.2 
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0         0.1        0.2        0.3  0.4 
AP(%) 

Figure  2. 10  Maximum  achievable  aspect  ratio  versus  pressurization  for  a 
Bond  number  of  8.8.  The  pressurization  is  calculated  as  a  percentage  of  the 
ambient  pressure,  Pambient  =  l-0x  105  Nm-2. 


2.6.4    Solutions  for  Aspect  Ratio  of  Unity 

Calculations  were  carried  out  for  aspect  ratio  of  unity,  Bond  number  of  8.8  and  six 
values  of  pressurization  in  the  range,  0.2  <  AP  <  0.3.  Figures  2.11  and  2.12  show  the 
results  of  the  calculations. 

Figure  2.11  shows  the  profiles  and  the  corresponding  free  energy  curves  that  are 
obtained  when  contact  angles  of  30  deg  and  60  deg  are  imposed.  The  free  energy  curves  also 
indicate  the  range  of  contact  angles  for  which  the  Young-Laplace  equation  can  be  integrated. 
It  is  useful  to  draw  attention  to  the  case  of  (j)0  =  30  deg,  AP  =  0.2295%.  It  may  be  observed 
that  two  profiles  are  obtained,  marked  in  Fig.  2. 1  l(i)  with  the  solid  curves,  one  of  them  also 
marked  with  'o'.  The  corresponding  free  energy  curve  in  Fig.  2. 11  (hi)  is  marked  with  '□' 
and  shows  that  two  solutions  can  be  obtained,  one  of  which  satisfies  a  local  minimum.  In 
Fig.  2.11(i),  the  unmarked  solid  curve  satisfies  a  local  minimum  indicating  static  stability. 
For  (|>0  =  60  deg,  none  of  the  obtained  profiles  satisfy  a  local  minimum. 
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Figure  2.11  Calculations  for  Bond  number  of  8.8,  unity  aspect  ratio  and 
imposed  contact  angle  condition. 

(i)  and  (iii)  Profiles  meeting  a  desired  contact  condition  of  60  deg  and  30  deg 
respectively,  (ii)  and  (iv)  Free  energy  curves  for  imposed  contact  angles  of  60  deg 
and  30  deg  respectively,  for  the  various  values  of  the  pressurization  considered. 

The  curves  in  (ii)  and  (iv)  use  the  same  convention  to  indicate  the  pressurization 
value,  (ii)  and  (iv)  also  indicate  the  range  of  values  of  (j)c  that  are  obtainable  by 
integrating  the  Young-Laplace  equation.  It  may  be  noted  that  for  <f)0  =  30  deg,  and 
AP  =  0.2295%,  there  are  two  profiles  satisfying  the  given  contact  condition, 
indicated  by  the  solid  curves.  The  profile  marked  with  o  does  not  minimize  the 
free  energy  whereas  the  unmarked  one  does.  None  of  the  other  profiles 
corresponds  to  a  local  minimum  indicating  that  although  they  mathematically 
satisfy  the  Young-Laplace  equation,  they  are  not  statically  stable. 
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Figure  2. 12  shows  profiles  that  (i)  meet  a  fixed  radius,  rc,  condition  for  three  values  of 
rc/rb  =  1.0,  0.75  and  0.5,  and  (ii)  meet  a  fixed  volume  constraint  of  V  =  2  and  3 
(non-dimensionalized).  The  corresponding  rc  versus  §c  and  V  versus  <\>c  show  the  overall 
range  over  which  the  Young-Laplace  equation  can  be  integrated.  It  may  also  be  noted  that 
the  overall  solution  behavior  does  not  show  any  abrupt  changes  as  the  pressurization,  AP, 
is  varied  over  the  specified  range.  Hence,  the  curves  shown  in  Fig.  2.12(iv)  help  in 
determining  the  existence  of  solutions  that  can  simultaneously  satisfy  a  given  contact  angle 
condition  and  a  fixed  volume.  The  free  energy  curves,  such  as  those  shown  in  Fig.  2. 12(iii), 
show  whether  the  calculated  profiles  shapes  are  statically  stable. 

2.7  Concluding  Remarks 

In  this  chapter,  a  formulation  has  been  presented  for  fluid  flow  and  heat  transfer  for 
thermally  induced  buoyancy  driven  flow  and  thermocapillary  convection.  This  chapter  also 
addresses  the  problem  of  integrating  the  Young-Laplace  equation  to  determine  the  free 
surface  shapes  formed  during  float  zone  growth.  Questions  regarding  the  existence  and 
uniqueness  of  the  solutions  to  the  static  Young-Laplace  equation  are  addressed  in  the  context 
of  a  free  energy  minimization  criterion.  The  method  is  validated  against  available 
experimental  data.  A  variety  of  calculations  have  been  carried  out  to  determine  free  surface 
shapes  over  a  range  of  Bond  numbers,  aspect  ratios,  internal  pressurization  and  different 
boundary  conditions.  These  calculations  form  the  building  block  for  determining  the 
meniscus  in  float  zone  calculations.  The  governing  equations  and  the  solution  methodology 
discussed  in  this  chapter  have  been  used  for  all  the  convective  heat  transfer  computations  in 
this  work. 
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Figure  2.12     Calculations  for  Bond  number  of  8.8,  unity  aspect  ratio  and 
either  (a)  imposed  radius  at  the  top  or  (b)  imposed  volume  of  the  melt. 
Profiles  for  (i)  fixed  volume  and  (iii)  fixed  rc  (for  three  different  rc). 

(ii)  volume  versus  (j)c  curves  and 
(iv)  rc/rb  versus  (j)c. 

The  curves  in  (ii)  and  (iv)  use  the  same  convention  to  indicate  the  pressurization 
value,  (ii)  and  (iv)  also  indicate  the  range  of  values  of  (}>c  that  are  obtainable  by 
integrating  the  Young-Laplace  equation. 


CHAPTER  3 

C0NVEC1TVE  HEAT  TRANSFER  COMPUTATIONS 
3.1  Background 

In  this  chapter,  the  implementation  of  various  convection  schemes  is  discussed  using  the 
benchmark  solution  for  natural  convection  in  a  cavity  as  the  vehicle  for  testing  and 
development.  The  developed  techniques  are  then  applied  to  compute  transient  natural 
convection  in  an  enclosure  and  thermocapillary  convection  in  a  liquid  bridge  of  molten  NiAl. 
The  validity  of  the  static  free  surface  assumption  of  the  previous  chapter  is  also  assessed  by 
an  a  posteriori  estimation  based  on  the  thermocapillary  convection  in  the  liquid  bridge. 
3.2  Convection  Treatment  and  Accuracy  Assessment 

Treatment  of  the  convection  terms  is  one  of  the  most  fundamental  issues  in 
Navier-Stokes  computations.  It  is  well  established  that  the  first  order  upwind  scheme  is  the 
most  stable  one  and  does  not  produce  numerical  oscillations,  at  the  cost  of  smearing  out  the 
sharp  gradients  which  are  usually  present  in  a  convection  dominated  flow.  To  overcome  the 
problem  of  excessive  numerical  diffusion,  numerous  treatments  have  been  proposed  which 
are  nominally  second  order  accurate,  based  on  a  Taylor  series  expansion,  but  cannot 
guarantee  physical  realizability.  Such  schemes  cannot  be  used  with  confidence  for  some 
scalar  variables  such  as  temperature,  turbulent  kinetic  energy  or  rate  of  turbulent  energy 
dissipation,  which  should  be  non-negative.  Hence,  it  is  not  uncommon  to  adopt  the  practice 
of  treating  the  convection  terms  in  the  momentum  equations  with  higher  order  schemes, 
while  employing  a  lower  order  scheme  for  transport  equations  of  conserved  scalar  quantities. 
One  such  low  order  scheme  is  the  hybrid  scheme  which  is  a  combination  of  the  first  order 
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upwind  scheme  and  the  second  order  central  difference  scheme.  The  impact  of  such  a  mixed 
discretization  on  the  solution  accuracy,  however,  needs  to  be  assessed. 

An  appropriate  test  vehicle  for  such  an  assessment  is  the  natural  convection  of  air  in  a 
square  cavity  with  differentially  heated  side  walls.  Accurate  benchmark  solutions  for  this 
problem  exist  in  literature  (de  Vahl  Davis  1983,  Hortmann  et.  al.  1990)  and  are  available  for 
comparison  with  the  present  calculations.  The  Rayleigh  number  for  the  calculations  is  106 
and  the  Prandtl  number  for  air  is  0.7 1 .  The  top  and  bottom  walls  are  insulated.  The  right  and 
left  walls  are  maintained  at  constant  temperatures  of  Tq  and  T\  respectively;  in 
nondimensional  terms  they  can  be  normalized  to  1 .0  and  0.0  respectively.  The  governing 
equations  are  the  two  dimensional  steady  Navier-Stokes  and  energy  equations  (2.9).  The 
Boussinesq  approximation  has  been  employed  wherein  the  density  variations  are  neglected 
in  all  the  terms  of  the  governing  equations  except  the  source  term  for  the  v-momentum 
equation.  Figure  3.13  shows  a  schematic  of  the  problem  domain. 
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y*  f£  =  0 

dy 


T,  =  1 


Tavg  =  0.5 


T0  =  0 
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Figure  3.13     Natural  convection  of  air  in  a  differentially  heated  square  cavity. 
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An  important  feature  of  the  solution  is  that  it  satisfies  the  following  centrosymmetry 
conditions. 

u(x,y)  =  -  u(l  -  x,  1  -  y) 
v(x,y)  =  -  v(l  -  x,  1  -  y) 

T(x,y)  =  -(T(l  -x,l  -y)-Tavg)  ^ 
il)(x,y)  =  1|>(1  -  x,  1  -  y) 

where  ip(x,  y)  is  the  stream  function  and  Tavg  =  0.5  (T,  +  T0).  While  this  feature  of  the 
solution  was  not  directly  utilized  in  the  computational  procedure,  it  was  used  as  an  indicator 
of  convergence  of  the  solutions  to  machine  accuracy. 

The  solution  technique  was  the  SIMPLE  algorithm  described  in  the  previous  section. 
The  central  difference  method  was  used  to  evaluate  the  diffusion  fluxes.  The  following 
convection  schemes  employing  different  profile  assumptions  were  tested  on  the  convection 
terms: 

(i)  Central  differences  for  both  momentum  and  energy  equations, 

(ii)  Second  order  upwind  for  both  momentum  and  energy  equations,  and 

(iii)  A  mixed  treatment  consisting  of  the  second  order  upwind  scheme  for  the  momentum 
equations  and  the  hybrid  scheme  for  the  energy  equation. 

It  may  be  noted  that  while  the  implementation  of  the  central  differencing  scheme  is  well 
established  due  to  its  simplicity,  that  of  the  second  order  upwind  is  not.  The  central  difference 
scheme  computes  the  convective  flux  as  a  linear  interpolation  of  the  values  in  two 
neighboring  control  volumes.  However,  the  second  order  upwind  scheme  has  been  treated 
in  a  variety  of  ways  resulting  in  inconsistent  performance  in  literature  (Shyy  et.  al.  1992, 
Vanka  1987).  The  scheme  adopted  in  this  study  is  that  proposed  by  Shyy  et.  al.  (1993)  where 
the  convective  flux  is  evaluated  by  a  linear  extrapolation  of  the  dependent  variable  from  two 
nearest  upwind  neighbors.  This  results  in  a  conservative  treatment  and  yields  consistent 
performance  with  respect  to  grid  refinement  for  two  dimensional  driven  cavity  flows  (Shyy 
et.  al.  1 993).  In  this  work,  the  same  implementation  is  tested  against  natural  convection  flow. 
The  present  investigations  and  the  results  of  Shyy  et.  al.  (1993),  together  give  a  more 
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comprehensive  account  of  the  performance  of  the  second  order  upwind  scheme  for  complex 
flow  calculations.  The  two  cases  represent  two  different  types  of  problems,  one  induced  by 
shear  stress  and  the  other  by  body  force. 

Two  benchmark  solutions  are  available  for  the  purpose  of  comparison.  One  is  due  to  de 
Vahl  Davis  (1983).  This  work  used  the  streamfunction  -  vorticity  formulation  of  the 
governing  equations.  Central  differencing  was  used  for  all  spatial  derivatives  and  the  method 
of  false  transients  using  first  order  forward  difference  in  time  was  used  to  solve  the  equations. 
The  finest  grid  used  was  an  uniform  81x81  mesh.  Richardson  extrapolation  was  then  used 
to  obtain  more  accurate  grid  independent  results  for  the  dependent  variables  and  the  derived 
quantities  such  as  the  Nusselt  number.  Another  benchmark  solution  is  that  due  to  Hortmann 
et.  al.  (1990).  This  work  solved  the  Navier-Stokes  equations  using  a  SIMPLE-like  solution 
algorithm  for  a  finite  volume  formulation  on  a  non-staggered  grid.  Central  differencing  was 
used  to  discretize  all  convection  terms.  A  multigrid  procedure  was  used  to  solve  the  resulting 
systems  of  linearized  equations.  Computations  were  carried  out  for  grid  sizes  up  to  321  X 
321  on  an  uniform  grid  and  up  to  641  X  641  on  a  non-uniform  grid.  The  results  are 
considered  accurate  to  within  0.01%. 

In  the  current  study,  computations  have  been  carried  out  on  uniform  grids  of  up  to 
81x81  using  a  variety  of  approximations  for  the  convection  terms.  Comparisons  have  been 
carried  out  with  the  solution  of  Hortmann  et.  al.  (1990)  for  all  quantities  considered  except 
the  streamfunction.  Since  the  streamfunction  data  is  not  available  in  Hortmann  et.  al.  (1990), 
the  required  data  has  been  obtained  from  de  Vahl  Davis  (1983). 

Computation  of  the  local  heat  transfer  rate  utilizes  the  staggered  grid  arrangement  for 
a  consistent  evaluation  of  the  heat  flux.  The  expression  for  the  convective  part  of  the  heat 
flux  is  consistent  with  the  convection  scheme  used  in  the  computations.  This  ensures  that  the 
total  heat  flux  crossing  any  vertical  line  is  a  constant  throughout  the  cavity. 
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3.2.1    Results  and  Discussion 

Figure  3.14  shows  the  contour  plots  of  the  streamfunction,  isotherms  and  heat  flux.  The 
streamfunction  plots  shows  the  vertical  boundary  layers  developing  along  the  isothermal 
walls  and  the  horizontal  boundary  layers  along  the  adiabatic  walls.  Centrosymmetry  of  the 
solution  can  be  observed.  Also,  three  secondary  vortices  rotating  in  the  clockwise  direction 
can  be  seen.  Two  of  these  are  displaced  towards  the  walls  in  the  streamwise  direction  and 
the  third  is  a  central  vortex  that  separates  these  two.  Distortion  of  the  isotherms  in  the 
streamwise  direction  can  be  observed  in  these  regions.  It  should  be  noted  that  the  heat  flux 
function  cannot  satisfy  the  centrosymmetry  condition. 

Comparisons  of  the  obtained  solutions  with  the  benchmark  solutions  have  been  made. 
The  following  quantities  have  been  used  for  comparison  and  error  analysis  in  this  study. 

(i)  Streamfunction  at  the  midpoint  of  the  cavity, 

(ii)  Average  Nusselt  number, 

(iii)  Maximum  u-velocity  on  the  vertical  mid-plane  of  the  cavity  (umax),  and 

(iv)  Maximum  v-velocity  on  the  horizontal  mid-plane  of  the  cavity  (vmax). 

Figure  3.15  shows  the  percentage  relative  errors  versus  the  grid  spacing,  h,  for  all  solutions 
due  to  the  different  convection  schemes  used.  It  can  be  observed  that,  for  the  derived 
quantities  —  the  Nusselt  number  and  the  streamfunction  —  the  mixed  treatment  is  the  most 
accurate.  Overall,  the  three  schemes  exhibit  comparable  performance  with  no  clearly 
superior  one  emerging.  The  errors  are  not  monotonic  with  decreasing  grid  spacing. 
However,  the  locations  of  the  quantities  umax  and  vmax  tend  to  vary  with  change  in  the  grid 
spacing  indicating  that  they  may  not  be  accurate  indicators  of  grid  independent  solutions. 
Nevertheless,  Fig.  3.15  clearly  shows  that  a  mixed  treatment  utilizing  the  second  order 
upwind  scheme  for  the  momentum  equations  and  the  hybrid  scheme  for  the  energy  equation 
yields  comparable  performance  with  the  second  order  central  difference  scheme. 

The  hybrid  scheme  automatically  switches  from  a  central  difference  to  a  first  order 
upwind  scheme  at  a  cell  Peclet  number  of  2.  The  cell  Peclet  numbers  are  defined  as, 
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(c)  Heat  function,  min:  -9.0  max:  9.0 
AH  =  0.9.  Contour  values  greater  than 
or  equal  to  0.0  have  been  thickened. 

Figure  3.14  Natural  convection  in  a  square  cavity.  Solutions  for  Ra 
=  106  and  Pr  =  0.71.  (a)  Streamf unction  (b)  Isotherms  (c)  Heat  flux 
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Figure  3.15  Percentage  error  in  solution  versus  grid  spacing,  h. 
(a)  umax  (b)  vmax  (c)  Streamfunction  (d)  Nusselt  number 

Key: 

  Central  difference 

  Mixed  treatment 

  Second  order  upwind 
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Peu  =  QCpuAy/k 
Pev  =  gCpvAx/k 

and  their  contours  are  shown  in  Fig.  3.16.  It  has  also  be  observed  that  if  the  cell  Peclet 
number  is  greater  than  2  and  the  solution  gradients  are  sharp,  then  the  central  difference 
solutions  may  be  wiggle-free  yet  inaccurate  (Shyy  et.  al.  1992).  A  good  measure  of  the 
solution  gradients  are  the  cell  Peclet  numbers  based  on  Au  and  Av  rather  than  u  and  v. 
Contours  of  the  cell  Peclet  number  based  on  Au  and  Av  are  also  shown  in  Fig.  3.16.  It  can 
be  observed  that,  on  the  81  X  81  grid,  Peu  is  less  than  1.6  over  the  entire  domain  indicating 
that  the  central  difference  scheme  is  used  for  the  u-momentum  equation  whereas  Pev 
exceeds  2  in  the  vertical  boundary  layers  close  to  the  isothermal  walls  indicating  that  the  first 
order  upwind  scheme  is  used  in  the  v  momentum  equations  inside  the  vertical  boundary 
layers.  It  should  be  noted  that  on  the  coarser  grids,  the  hybrid  scheme  will  switch  to  the  first 
order  upwind  scheme  in  the  boundary  layer  regions  due  to  the  higher  values  of  Ax  and  Ay. 
Nevertheless,  Fig.  3. 16  indicates  that  for  all  grid  spacings  tested,  the  use  of  the  lower  order 
convection  treatment  in  the  energy  equation  does  not  degrade  the  solution  accuracy. 
Thus,  it  is  clearly  established  that 

(i)  the  second  order  upwind  scheme  yields  satisfactory  solutions,  and 

(ii)  an  apparently  lower  order  combination  of  schemes  consisting  of  second  order  upwind 
for  the  momentum  equations  and  the  hybrid  scheme  for  the  energy  equation  can  produce 
satisfactory  results  for  this  problem. 

3.3  Transient  Natural  Convection  in  an  Enclosure 
In  the  following,  the  present  numerical  method  will  be  further  validated  in  the  context 
of  a  transient  natural  convection  flow  relevant  to  electronics  cooling.  In  the  literature,  two 
generic  problem  types  have  been  extensively  investigated;  the  enclosure  and  the  open  cavity. 
Flows  within  enclosures  are  bounded  on  all  sides  by  walls,  whereas  open  cavities,  such  as 
vertical  parallel  plate  channels  or  tubes,  have  their  ends  open  to  the  ambient  flow.  An 
extensive  discussion  of  important  heat  transfer  results  pertaining  to  such  configurations  has 
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(a)  Cell  Peclet  numbers  based  on  u.  (b)  Cell  Peclet  numbers  based  on  v.  min: 
min:-1.6(-)max:  1.6  (+)  -2.8  (-)  max:  2.8  (+)  APev  =  0.7. 

APeu  =  0.4.  Contour  values  less  than  -2  or  greater 

than  2  have  been  thickened. 


(c)  Cell  Peclet  numbers  based  on  Au. 
min:  -0.22  max:  0.18  APeAu  =  0.04. 
Contour  values  greater  than  or  equal  to 
zero  have  been  thickened. 


(d)  Cell  Peclet  numbers  based  on  Av.  min: 
-0.18  max:  0.22  APeAv  =  0.04.  Contour 
values  greater  than  or  equal  to  zero  have 
been  thickened. 


Figure  3.16     Contours  of  cell  Peclet  number. 

(a)  based  on  u  (b)  based  on  v  (c)  based  on  Au  (d)  based  on  Av. 
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been  presented  by  Raithby  and  Holland  ( 1 985),  Gebhart  et  al  ( 1 988)  and  Peterson  and  Ortega 
(1990). 

The  first  experimental  investigation  of  natural  convection  of  air  in  isothermal, 
parallel-plate  channels  was  carried  out  by  Elenbaas  (1942).  By  analyzing  a  simplified  set 
of  equations,  and  fitting  experimental  results,  a  correlation  for  the  Nusselt  number  was 

obtained  against  Ra*  where  Ra*  =  Ra  d/L  is  the  Elenbaas  Rayleigh  number,  Ra  is  the 
Rayleigh  number  based  on  channel  half  width,  d  is  the  half  width  of  the  channel  and  L  is  its 

height.  The  experiments  spanned  the  range  0. 1  <  Ra*  <  105.  Subsequent  investigations 
by  Sparrow  and  Bahrami  (1980)  showed  inaccuracies  in  the  data  of  Elenbaas  at  low  values 
of  Ra*  and  they  observed  that  the  corrections  employed  by  Elenbaas  for  extraneous  heat 
losses,  like  radiation,  exceeded  the  actual  natural  convection  heat  transfer;  they  and  Martin 
et  al  (1991)  have  improved  Elenbaas 's  correlation  for  the  isothermal  parallel  plate  channel. 

Isothermal,  asymmetrically  heated,  vertical  parallel-plate  channels  have  been 
investigated  by  Sparrow  et  al  (1984).  In  this  study,  one  wall  was  isothermal  while  the  other 
was  unheated.  The  experiments,  carried  out  in  water,  with  Pr  =  5,  revealed  a  pocket  of 
recirculating  flow  adjacent  to  the  unheated  wall,  close  to  the  downstream  end  of  the  channel. 
In  spite  of  this  recirculating  zone,  the  Nusselt  number  correlated  well  with  the  parameter 
(2d/L)  Ra  where  2d  is  the  inter-wall  spacing.  They  also  complemented  the  experimental 
results  with  numerical  solutions  based  on  parabolic  equations,  derived  by  neglecting  the 
streamwise  second  derivatives.  Thus,  they  concluded  that  the  Nusselt  number  is  relatively 
insensitive  to  the  backflow.  Sparrow  and  Azevedo  (1985)  followed  up  this  study  by 
investigations  over  a  greater  range  of  inter-wall  spacing.  They  were  able  to  obtain  a 
correlation  for  the  Nusselt  number  against  (2d/L)  Ra  over  a  wide  range  spanning  the  fully 
developed  limit  and  the  flat  plate  limit.  They  found  that  the  flat  plate  limit  does  not  form 
an  upper  bound  on  the  heat  transfer  results.  They  also  concluded  that  the  Nusselt  number 
is  insensitive  to  3-D  effects. 
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Experiments  have  also  been  performed  for  iso-heat-flux  plates  by  Wirth  and  Stutzman 
(1982)  for  the  symmetrically  heated  case  and  by  Webb  and  Hill  (1989)  for  the 
asymmetrically  heated  plates.  Correlations  were  presented  for  the  average  Nusselt  number, 
local  Nusselt  number  and  the  maximum  temperature  on  the  plates,  as  a  function  of  the 
Rayleigh  number.  The  results  were  in  reasonable  agreement  with  the  previous  numerical 
work  of  Aung  (1972),  for  the  limit  of  fully  developed  flow,  and  Aung  et.  al.  (1972)  for 
developing  flow. 

Theoretical  investigation  of  free  convection  between  parallel-plates  was  first 
undertaken  by  Bodoia  and  Osterle  (1962).  They  solved  the  finite  difference  form  of  the 
parabolic  equations,  and  obtained  good  agreement  with  the  data  of  Elenbaas  except  at  low 
values  of  the  Grashof  number,  which  they  attributed  to  side  leakage  effects.  However,  they 
confined  their  analysis  to  the  slot  region,  neglecting  the  effects  of  inflow  and  outflow 
boundaries. 

In  order  to  resolve  the  flow  profile  at  the  inlet,  Kettleborough  (1972)  solved  the  laminar, 
2-D  Boussinesq  equations  for  air  with  Pr  =  0.7,  Gr  =  12.5  and  1250,  and  on  a  domain  with 
an  aspect  ratio,  PR  of  5  that  included  an  upstream  entrance  region.  He  however,  imposed 
fully  developed  conditions  at  the  channel  exit,  which  is  a  restrictive  assumption.  A  similar 
analysis  was  carried  out  by  Nakamura  et.  al.  (1982)  with  the  difference  that  the  wall 
boundary  condition  for  the  streamfunction  was  obtained  by  consideration  of  the  pressure 
drop  along  the  channel  wall.  Their  results  were  significantly  different  from  those  of 
Kettleborough,  especially  at  higher  Grashof  number.  Also  they  could  not  capture  the  reverse 
flow  region  predicted  by  Kettleborough. 

Chang  and  Lin  (1989)  performed  a  transient  analysis  which  ultimately  ended  in  steady 
state.  Calculations  were  carried  out  for  Pr  =  0.7,  103  <  Ra  <  106,  and  AR  =  5,  10.  The  inlet 
and  exit  domains  were  considered  in  the  analysis.  However,  they  replaced  the  180°  turn  at 
the  slot  inlet  and  exit  with  an  artificial  90°  adiabatic  corner.  They  presented  the  time 
evolution  of  the  flow  and  thermal  fields  and  also  a  correlation  for  the  Nusselt  number  as  a 
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function  of  the  Rayleigh  number  and  the  channel  aspect  ratio,  for  air.  At  the  free  surfaces 
they  prescribed  the  following  boundary  conditions;  normal  gradient  of  all  velocities  were 
zero  at  the  free  surfaces,  normal  gradient  of  temperature  was  zero  if,  locally,  the  fluid  flowed 
out  of  the  domain,  otherwise,  the  temperature  was  set  equal  to  the  ambient  temperature. 

Naylor  et  al  (1991)  have  solved  the  steady  2-D  Navier-Stokes  equations  with 
consideration  of  the  inlet  domain.  They  also  replaced  the  180°  turn  at  the  inlet  and  outlet 
with  an  artificial  90°  adiabatic  corner.  The  inflow  boundary  was  a  semicircle  and  the 
boundary  conditions  were  obtained  from  considerations  of  Jeffrey-Hamel  flow  which  is  a 
similarity  solution  of  isothermal  flow  due  to  a  sink  at  the  point  of  intersection  of  two  walls. 
Fully  developed  flow  conditions  were  imposed  at  the  slot  exit.  They  also,  for  comparison 
purposes,  solved  the  parabolic  equations.  Their  solutions  are  in  closer  qualitative  agreement 
with  those  of  Nakamura  et  al  (1982)  than  with  Kettleborough  (1972),  but  the  overall 
quantitative  agreement  with  either  is  quite  poor.  Reverse  flow  predicted  by  Kettleborough 
was  not  obtained. 

Shyy  et  al  (1992)  carried  out  calculations  modelling  both  the  inflow  and  the  outflow 
domains.  The  inflow  domain  had  an  upstream  boundary  condition  amounting  to  that 
generated  by  a  virtual  line  mass  source  located  at  the  centerline  of  the  inlet,  and  whose  side 
boundaries  were  modelled  by  side  entrainment  boundary  conditions.  For  the  outflow  region, 
the  downstream  boundary  corresponded  to  the  plume  solution  generated  by  a  virtual  line  heat 
source  located  at  the  inlet  centerline.  The  source/sink  strengths  were  determined  by 
matching  the  total  convective  energy  transported  across  each  boundary  with  the  thermal 
energy  conducted  from  the  plates  to  the  fluid.  For  the  side  entrainment  boundaries,  the 
normal  derivatives  of  the  velocity  components  were  set  to  zero  and  the  normal  derivative  of 
the  temperature  was  set  to  zero,  if  the  fluid  flowed  out  of  the  domain;  otherwise,  the 
temperature  was  set  equal  to  the  ambient  temperature.  Calculations  were  carried  out  for  Gr 
=  103  ,  104  ,  105  and  for  AR  =  1  ,  2.  For  channel  aspect  ratio  of  unity,  and  at  the  Grashof 
numbers  studied,  a  recirculation  zone  was  observed  in  the  slot. 
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Ramanathan  and  Kumar  ( 1 99 1 )  tried  to  circumvent  the  problem  of  artificial  outflow  and 
inflow  boundary  conditions  by  putting  the  slot  within  a  large,  isothermal  enclosure. 
Iso-heat-flux  boundary  conditions  were  applied  to  the  side  of  the  heated  plates  facing  the 
interior  of  the  channel.  The  sides  facing  the  enclosure  were  insulated.  The  problem  is 
difficult  to  solve  parametrically;  apart  from  the  usual  Ra,  Pr,  and  channel  AR ,  there  are  five 
more  geometric  parameters.  They  found  that  their  results  were  sensitive  to  the  relative 
location  of  the  slot  in  the  vertical  direction.  In  order  to  extrapolate  their  results  to  the  limit 
of  infinitely  large  enclosures,  the  temperatures  within  the  slot  were  offset  by  the  entry 
temperature  at  the  inlet  centerline.  The  steady  2-D  Navier-Stokes  equations,  incorporating 
the  Boussinesq  approximation,  were  employed  for  the  computations.  They  obtained  good 
agreement  with  the  experimental  results  of  Wirth  and  Stutzman  (1982)  for  Ra  =  18.8,  191, 
and  2414.  Correlations  were  presented  for  maximum  temperature  on  the  plates  and  for  the 
average  Nusselt  number  versus  the  Rayleigh  number  and  channel  aspect  ratio,  for  air.  They 
also  found  that  for  fluids  with  Pr  >  0.7,  the  hydrodynamic  development  length  is  smaller  than 
the  thermal  development  length  up  to  Ra  =  400  and  for  Pr  <  0.1,  the  opposite  holds.  They 
concluded  that  vertical  diffusion  of  energy  cannot  be  neglected,  especially  for  low  Pr  and 
short  /R  channels. 

Thus  it  is  seen  that  various  investigators,  in  the  attempt  to  compute  natural  convection 
heat  transfer  between  two  vertical  plates  have  formulated  various  boundary  conditions  to 
deal  with  the  artificial  boundaries  of  a  finite  computational  domain.  The  results  are  also  as 
varied.  The  heat  transfer  results  seem  to  correlate  well  with  parameters  obtained  by 
neglecting  streamwise  diffusion,  yet,  in  terms  of  the  velocity  field,  recirculating  flows  have 
been  found  within  the  slot,  which  is  incompatible  with  the  boundary  layer  assumption. 
Attempts  have  been  made  to  solve  the  full  elliptic  equations  with  the  fully  developed 
boundary  conditions  at  the  channel  exit,  which  is  restrictive  and  has  been  found 
unsatisfactory  by  other  investigators. 
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Despite  the  fact  that  most  investigators  have  dealt  with  the  flow  between  two  parallel 
plates  in  an  open  domain,  many  practical  applications,  especially  in  the  electronics  industry, 
have  geometrical  configurations  similar  to  the  one  in  this  study,  as  shown  in  Fig  3.17.  For 
example,  as  discussed  by  Gingrich  et  al  (1992),  an  electronics  cooling  cabinet  often  adopts 
such  a  generic  design  where  hot  plates  are  placed  in  the  middle  of  the  enclosing  cabinet.  Such 
geometries,  however,  involve  several  independent  geometrical  parameters  and  it  may  not  be 
feasible  to  conduct  parametric  investigations  over  the  entire  range  of  possible  values, 
especially  if  the  transient  effects  are  of  interest.  Recently,  interest  has  arisen  in  the  use  of 
non-Newtonian  fluids  in  electronics  cooling  applications  (Gingrich  et  al  1992),  which, 
however,  is  outside  the  scope  of  the  current  study. 

A  discussion  of  the  present  contribution  relative  to  the  previous  efforts  is  in  order.  In 
Shyy  et.  al.  (1992),  not  only  were  the  computations  conducted  for  an  open  channel,  but  more 
importantly,  there  was  a  major  issue  in  terms  of  prescribing  the  boundary  conditions  for  that 
configuration.  As  discussed  in  that  paper,  there  is  a  fundamental  difficulty  in  determining 
the  overall  incoming  and  outgoing  mass  fluxes  and  their  distributions,  for  an  open  channel 
flow,  unless  some  simplifications  are  made  in  the  problem  formulation.  That  is  the  rationale 
behind  deriving  the  line-source  and  the  plume  solutions  for  the  inlet  and  outlet  conditions 
respectively.  A  detailed  account  of  this  issue  can  be  found  in  Shyy  et.  al.  (1992).  There  is 
a  need  to  understand  the  way  of  extending  the  theoretical  framework  from  a  confined 
enclosure  to  that  of  an  open  domain,  and  the  correct  approach  has  yet  to  be  identified.  The 
configuration  chosen  here  represents  a  more  realistic  case  in  the  area  of  electronics  cooling. 
Furthermore,  the  present  work  is  a  simulation  of  transient  heat  transfer,  whereas  the  work 
of  Shyy  et.  al.  (1992)  is  a  steady  state  formulation. 

The  current  study  consists  of  two  isothermal  vertical  parallel  plates  placed  within  a 
rectangular  enclosure.  The  top  and  bottom  walls  of  the  enclosure  are  insulated,  whereas  the 
side  walls  are  isothermal,  as  shown  in  Fig  3.17.  The  geometrical  parameters  are  held 
constant    and    have    the   values    Xmax/2d  =  4,  Ymax/2d  =  3,  L,/2d  =  1.5  and 
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Figure  3.17  Schematic  of  present  problem  with  a  heated  vertical  channel  inside 
a  cold  rectangular  enclosure.  The  hatched  lines  indicate  adiabatic  walls.  The 
domain  is  symmetric  about  the  vertical  centerline.  The  side  walls  of  the  enclosure 
and  the  channel  walls  are  isothermal.  The  letters  A,B  indicate  the  vertical  extent 
of  the  enclosure  wall.  Similarly,  the  letters  C,D  designate  the  vertical  extent  of  the 
channel  wall. 

L2/2d  =  1.  Initially  there  is  no  convection  and  uniform  temperature  prevails.  At  time  t  = 
0,  the  plates  experience  a  step  change  in  the  temperature.  A  Grashof  number  of  105,  typical 
of  electronic  devices,  is  considered  in  our  work.  The  present  study  attempts  to  accurately 
capture  the  detailed  flow  and  thermal  structures  of  the  transient,  two-dimensional  buoyancy 
induced  flowfield  and  associated  heat  transfer  results.  Thus,  care  has  been  taken  to  ensure 
adequate  numerical  accuracy  in  both  time  and  space.  Also  reasonable  for  the  Navier-Stokes 
flow  fields,  is  the  Boussinesq  approximation,  which  has  been  employed  here.  The  main 
motivation  of  the  work  is  to  utilize  a  set  of  accurate  transient  simulations  to  facilitate  a 
detailed  scrutiny  of  the  development  of  the  flow  and  heat  transfer  characteristics,  yielding 
improved  insights  in  important  practical  problems. 
3.3.1    Formulation  and  Solution  Procedure 

The  non-dimensionalized,  unsteady,  two-dimensional,  Navier-Stokes  equations 
incorporating  the  Boussinesq  approximation  can  be  written  as: 
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The  dimensionless  variables  are  defined  by  Eq.  (2.5),  where  the  reference  scales  are  d,  the 
half  channel  width,  the  velocity  scale  Ur  =  a/d,  a,  the  thermal  diffusivity,  and  the  time 

scale  tr  =  d2/a.  The  Boussinesq  approximation  involves  the  assumption  that  density 
variations  are  small  enough  to  be  neglected  everywhere  except  in  the  buoyancy  source  term 
in  the  Y-momentum  equation.  All  other  fluid  properties  such  as  viscosity,  thermal 
conductivity  and  specific  heat  are  assumed  constant. 

The  solution  methodology  involved  the  use  of  the  second  order  central  differences  for 
discretization  of  all  spatial  derivatives,  including  the  convection  terms,  and  the  first  order, 
implicit,  backward  Euler  scheme  for  the  time  derivatives.  The  resulting  algebraic  equations 
were  iterated  till  a  desired  accuracy  was  achieved.  The  SIMPLE  algorithm,  Patankar  (1980), 
which  is  a  semi-implicit  pressure  correction  based  approach,  was  used  for  the  iteration 
procedure  at  each  time  step.  The  following  quantities  were  monitored  at  each  time  step:  (i) 
heat  flux  from  the  heated  plates,  (ii)  heat  flux  out  of  the  enclosure  wall,  and  (iii)  volume  flux 
through  the  channel.  These  quantities  were  nondimensionalized  as: 
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where  xp  =  Li  /  d  is  the  non-dimensional  x-location  of  the  heated  plate  under  consideration 
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where  x  =  0  defines  the  cold  enclosure  wall,  and 
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where     is  any  y-location  within  the  channel  region. 

The  temperature  gradients  at  the  walls  were  evaluated  by  one-sided  finite  differences 
at  the  walls.  The  integrals  were  evaluated  by  assuming  a  piece-wise  constant  distribution 
of  the  dependent  variables  within  each  control  volume,  which  is  consistent  with  the 
discretization  procedure.  The  heat  flux  distribution  on  the  isothermal  walls  can  be  described 
as  a  local  Nusselt  number: 
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where  xp  is  the  x-location  of  the  heated  plate,  and  x  =  0  defines  the  cold  enclosure  wall. 
Calculations  were  carried  out  for  a  Grashof  number,  Gr  =  105,  and  channel  aspect  ratio, 
AR  =  1.  Other  geometrical  parameters  were  also  held  constant. 

Initial  computations  were  carried  out  on  a  coarse  grid  of  60  X  60  control  volumes  and 
a  fine  grid  of  90  X  90  control  volumes,  with  a  time  step  (nondimensional)  At  =  10-2.  The 
grid  spacing  was  adjusted  by  a  cubic  polynomial,  with  the  grid  spacing  specified  at  the  two 
ends.  The  coarse  grid  had  30  control  volumes  between  the  enclosure  wall  and  the  channel 

wall,  in  the  x -direction,  with  AX/2d  =  10  ~  3  at  the  channel  wall  and  10~2  at  the  enclosure 
wall  and  30  control  volumes  between  the  channel  wall  and  the  channel  centerline.  In  the  y 
-  direction,  there  were  15  control  volumes  between  the  bottom  enclosure  wall  and  the 
channel  inlet,  20  control  volumes  in  the  channel  region  and  25  control  volumes  between  the 
channel  outlet  and  the  enclosure  top  wall.  The  grid  spacing  was  AY/2d  =  10~3  at  the  channel 
inlet  and  outlet  and  10~2  at  the  enclosure  wall.  The  corresponding  parameters  for  the  fine 
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grid  were  AX/2d  =  5  x  10-4  at  the  channel  wall  and  5  X  10-3  at  the  enclosure  wall,  with 
50  control  volumes  in  the  x  -direction  between  the  enclosure  side  walls  and  the  channel 
walls,  and  40  control  volumes  between  the  channel  wall  and  the  channel  centerline.  In  the 
y  -  direction,  the  parameters  were  AY/2d  =  5  X  10"4  at  the  channel  inlet  and  outlet,  and  5 
X  10-3  at  the  enclosure  top  and  bottom  walls.  There  were  20  control  volumes  in  the  region 
between  the  bottom  enclosure  wall  and  the  channel  inlet,  30  control  volumes  in  the  channel 
region,  and  40  control  volumes  in  the  region  between  the  channel  outlet  and  the  top  enclosure 
wall.  Symmetry  boundary  conditions: 

"  -  0  •  £  -  0  •  f  =  0  e»> 

were  applied  at  the  channel  centerline  and  the  calculations  were  carried  out  for  half  the 
domain,  thereby  achieving  good  spatial  resolution  along  with  economy  of  storage  and  CPU 
time.  At  each  time  step,  the  algebraic  equations  were  iterated  until  the  normalized  residuals 
obtained  by  summing,  over  the  whole  domain,  the  absolute  values  of  the  mass,  momentum 
and  energy  flux  imbalances  within  each  cell,  were  within  10-3. 
3.3.2    Results  and  Discussion 

Initially,  a  time  step,  At  =  10~2  was  chosen  and  the  calculations  advanced  up  to  300  time 
steps.  At  that  point,  the  calculations  converged  in  less  than  3  iterations  per  time  step.  At  that 
time  instant,  the  heat  balance  within  the  domain  was  within  0.6%  for  the  coarse  grid  and 
0.4%  for  the  fine  grid,  yielding  a  difference  of  less  than  0.2%  between  the  two  solutions.  The 
iterative  calculation  procedure  at  a  given  time  step  was  terminated  when  the  normalized 
residual,  for  each  dependent  variable,  was  less  than  10~3.  The  residual  for  any  given 
dependent  variable  is  computed  by  summing  over  the  absolute  values  of  the  relevant  flux 
imbalances  in  each  control  volume  over  the  entire  computational  domain.  The  residuals  are 
then  normalized  with  respect  to  characteristic  quantities  obtained  from  the  reference  scales 
defined  previously.  For  example,  the  mass  imbalance  in  the  domain  is  non-dimensionalized 
with  respect  to  the  total  mass  in  the  enclosure.  From  Fig.  3.18  it  can  be  observed  that  the 
above  comparisons  hold  true  over  the  entire  time  history.  It  was  clear  that  grid  independence 
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(c)  time  (d)  time 


Figure  3.18  Comparison  of  calculations  on  the  61  x  61  grid  versus  the  91  x  91 
grid.  At  =  10~2  for  both  cases. 

(a)  Number  of  iterations  for  convergence  at  each  time  step. 

(b)  Nusselt  number,  Nup,  based  on  heat  flux  through  channel  wall. 

(c)  Nusselt  number,  Nuc,  based  on  heat  flux  through  enclosure  wall. 

(d)  Nondimensional  volume  flux  through  channel,  Q. 

Key:   91  x  91  grid    61  x  61  grid. 


had  been  achieved.  However,  to  maintain  adequate  computational  accuracy,  the  fine  grid 
was  used  for  all  further  calculations. 

To  verify  the  time  step  independence  of  the  solution,  a  computation  was  conducted  for 
At  =  10 ~4  and  marched  up  to  104  time  steps.  Figure  3.19  summarizes  some  of  the  key 
findings  in  terms  of  the  monitored  quantities.  At  time  t  =  1.0,  the  differences  between  the 
two  solutions  are  within  0.2%  for  the  heat  flux  from  the  heated  plates,  0.3%  for  the  heat  flux 
out  of  the  enclosure  walls  and  0. 1  %  for  the  volume  flux  through  the  channel.  However,  at 
intermediate  time,  i.e.,  0.01  <  t  <  0.6,  there  are  significant  differences  between  the 
solutions.  For  the  Nusselt  number  based  on  heat  flux  from  the  channel  wall,  Nup,  the 
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Figure  3.19  Comparison  of  calculations  with  (a-c)  At  =  10~2  and  At  =  10"4  on  the 
91  x  91  grid  and  (d)  -  (f)  At  =  10"4  and  At  =  10~5  on  the  91  x  91  grid. 

(a)  Nusselt  number,  Nup,  based  on  heat  flux  through  channel  wall. 

(b)  Nusselt  number,  Nuc,  based  on  heat  flux  through  enclosure  wall. 

(c)  Nondimensional  volume  flux  through  channel,  Q 

(d)  Nusselt  number,  Nup,  based  on  heat  flux  through  channel  wall. 

(e)  Nusselt  number,  Nuc,  based  on  heat  flux  through  enclosure  wall. 

(f)  Nondimensional  volume  flux  through  channel,  Q 

(a)  -  (c)  Key:   At  =  10"4,  ■  At  =  lO"2 

(d)  -  (0  Key:   At  =  lO"5,  •  At  =  10"4 
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differences  are  negligible  throughout  the  time  history,  but  for  the  Nusselt  number  based  on 
heat  flux  out  of  the  enclosure  walls,  Nuc,  there  are  noticeable  quantitative  and  qualitative 
differences.  For  the  volume  flux,  Q,  the  solutions  corresponding  to  each  of  the  two  time  steps 
are  profoundly  different.  The  solution  corresponding  to  At  =  10-2  shows  a  highly  damped 
oscillatory  trend,  extending  up  to  t  ~  0.3,  but  does  not  show  net  mass  flow  reversal  through 
the  channel.  But  the  solution  corresponding  to  At  =  10"4  shows  oscillations  of  much  larger 
amplitude  showing  net  mass  flow  reversal  through  the  channel.  The  amplitude  reaches  a 
maximum  at  about  t  —  0. 1  and  the  oscillations  are  damped  to  their  steady  state  value  at  about 
t  ~  0.6.  The  peak-to-peak  average  also  displays  an  oscillatory  tendency.  The  solutions  also 
differ  in  phase.  As  an  explanation,  it  is  recalled  that  the  backward  Euler  is  first  order  in  time 
and  prone  to  numerical  damping.  As  a  check  on  the  solution  accuracy,  a  computation  was 
initiated  with  At  =  10~5  and  advanced  up  to  104  time  steps.  Figure  3.19  shows  no  discernible 
differences  in  the  solutions  corresponding  to  At  =  10-4  and  At  =  10~5,  indicating  that  the 
solution  is  independent  of  the  value  of  At. 

It  is  noted  that  the  solution  corresponding  to  At  =  10""4  was  marched  up  to  104  time  steps, 
i.e.,  t  =  1.0.  From  a  comparison  with  the  steady  state  solution  obtained  with  At  =  10~2,  it 
is  seen  that  steady  state  cannot  be  achieved  at  t  =  1 .0.  But  the  differences,  in  terms  of  the 
monitored  quantities  are  quite  small,  less  than  2.6%  for  Nup,  5%  for  Nuc  and  1%  for  Q.  The 
fine  time  step  calculations  are  expensive  in  terms  of  computer  time.  The  solution 
corresponding  to  At  =  10-4  and  104  time  steps  took  about  1.5  days  of  CPU  time  on  a 
DECstation  5000.  Since  the  steady  state  solution  can  be  obtained  by  using  At  =  10~2  or 
higher,  the  calculations  were  terminated  at  this  stage.  All  further  discussion  applies  to  the 
transient  solution  corresponding  to  At  =  10"4. 

Figure  3.20  shows  the  distribution  of  the  local  Nusselt  number  on  the  heated  channel 
walls.  The  solid  line  shows  the  Nusselt  number,  Nu+,  based  on  the  heat  flux  from  the  plates 
to  the  interior  of  the  channel  and  the  broken  line  the  Nusselt  number,  Nu-  based  on  the  heat 
flux  from  the  heated  plates  to  the  exterior  of  the  channel.  The  quantities  Nu+  and  Nu-  are 
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Figure  3.20  Time  evolution  of  Nusselt  number  distribution  on  the  inner  (Nu+) 
and  outer  (Nu-)  channel  wall  as  defined  in  Eq.  (3.8).  The  letters  C  and  D  designate 
the  channel  wall  at  2d  <  Y  <  4d,  X  =  1.5d ,  as  shown  in  Fig.  3.1. 

Key:   Nu+,  based  on  heat  flux  to  channel  interior 

 Nu-,  based  on  heat  flux  to  the  open  side. 
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defined  in  Eq.  (3.10).  These  are  established  early  in  time  and  do  not  vary  significantly 
thereafter.  Notably,  these  are  relatively  insensitive  to  the  oscillations  in  the  net  volume  flux 
through  the  channel,  Q,  as  will  be  discussed  later  in  connection  with  the  detailed  analysis  of 
the  flow  and  thermal  structures  in  the  slot. 

Figure  3.21  shows  the  local  Nusselt  number  distribution,  Nu  as  defined  in  Eq.  (3.10), 
on  the  enclosure  walls.  Initially,  there  is  a  spike-like  distribution  established  near  the  top, 
as  the  heated  fluid  reaches  that  location.  Gradually,  it  broadens  as  the  convection  field 
develops  and  increases  in  magnitude  and  at  large  time,  it  has  a  triangular  distribution 


Figure  3.21  Time  evolution  of  the  Nusselt  number  distribution  on  the 
enclosure  wall,  Nu  as  defined  in  Eq.  (3.8).  The  letters  A  and  B  refer  to  the 
vertical  extent  of  the  enclosure  wall  atO  <  Y  <  6d,  X  =  0  ,  as  shown  in 
Fig.  3.1 
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resulting  from  the  fact  that  the  convection  field  moves  towards  the  top  part  of  the  enclosure 
side-wall  and  moves  away  from  its  bottom  part. 

Figure  3.22  shows  the  streamlines,  isotherms  and  contours  of  motion  pressure  at  t  = 
0.01 , 0.05, 0. 1 , 0. 1 5, 0.2  and  1 .0.  Physically,  at  the  early  stages  of  the  transient  process,  heat 
transfer  from  the  heated  plates  to  the  fluid  amounts  to  a  conduction  process  (Gebhart  et  al. 
1988);  buoyancy  induced  motion  arises,  but  convective  effects  on  the  heat  transfer  process 
are  absent.  This  phenomenon  is  however,  present  only  for  a  short  time  during  the  startup 
of  the  transient.  At  any  location,  the  conduction  effect  ends  with  the  arrival  of  the  leading 
edge  effect.  Thereafter,  the  process  is  multi-dimensional  and  the  transient  process  continues 
towards  a  steady  state.  Thus  it  is  evident,  from  the  streamlines  and  the  isotherms,  that  time 
t  =  0.01  depicts  the  end  of  the  one  dimensional  conduction  process.  The  flow  adjacent  to 
the  channel  walls  consists  of  a  thin  layer  of  heated  fluid  moving  upward  and  has  a  structure 
qualitatively  similar  to  a  buoyancy  driven  boundary  layer  on  a  single  isothermal  vertical  flat 
plate  as  can  be  observed  from  the  v- velocity  profiles  across  the  middle  of  the  channel, 
depicted  in  Fig  3.23.  It  is  evident  from  Fig.  3.23  that  in  the  course  of  the  transient  evolution 
while  the  v-velocity  reverses  in  sign  in  the  regions  away  from  the  heated  plates,  adjacent  to 
the  plates  it  exhibits  only  modest  variations  in  magnitude  and  maintains  its  direction.  Thus 
the  flow  structure  near  the  channel  walls  is  qualitatively  insensitive  to  the  flow  in  the  rest 
of  the  channel.  This  particular  flow  feature,  established  early  in  the  transient  stages,  persists 
till  steady  state,  and  shields  the  wall  heat  transfer  process  from  the  effects  of  the  flow  in  the 
central  region  of  the  channel.  This  explains  the  relative  insensitivity  of  the  plate  Nusselt 
number  to  the  oscillating  volume  flux  within  the  slot.  At  t  =  0.05  large  vortical  structures 
can  be  seen  both  inside  and  outside  the  slot.  The  net  volume  flow  rate  in  the  channel,  Q,  is 
in  the  negative  y  direction.  The  isotherms  show  distorted  shapes  corresponding  to  the 
convection  field  in  the  vicinity  of  the  recirculating  flow  structures. 

At  t  =  0. 1  there  are  many  pockets  of  recirculating  flow  close  to  the  enclosure  walls  and 
near  its  top.  The  net  volume  flux,  Q,  is  directed  upward  and  attains  its  global  maximum 
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(a)  streamlines:  min:  -12.13  (between  plates) 
max:  21.15,  AW:  1.664 


(b)  isotherms:  min:  0  max:  0.997183 
AT:  0.0665 


(a)  Streamlines:  min:  -12.13  (between  plates) 
max:  21.15,  AW:  1.664 


(b)  isotherms:  min:  0  max:  0.997255 
AT:  0.0665 


high  motion  pressure 


(c)  motion  pressure:  AP:  421  (c)  motion  pressure:  AP:  2339 

Fig.  3.6(i)  t  =  0.01  Fig.  3.6(h)  t  =  0.05 

Figure  3.22  Contour  plots  of  the  solution. 
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(a)  Streamlines:  min:  -37.92  (-),  max:  90.25 
AW:  6.41 


(b)  isotherms:  min:  0  max:  0.997488 
AT:  0.0665 


(c)  motion  pressure:  AP:  1564 
Fig.  3.6(i)  t  =  0.1 

Figure  3.22  —  continued. 


(a)  Streamlines:  min:  -52.53  (between  plates) 
max:  99.1,  AW:  7.58 


(b)  isotherms:  min:  0  max:  0.99786 
AT:  0.0499 


(c)  motion  pressure:  AP:  2090 
Fig.  3.6(h)  t  =  0.15 
Contour  plots  of  the  solution 
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(c)  motion  pressure:  AP:  3229  (c)  motion  pressure:  AP:  6824 

Fig.  3.6(i)  t  =  0.2  Fig.  3.6(ii)  t  =  1.0 

Figure  3.22  —  continued.  Contour  plots  of  the  solution 
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value.  It  is  noted  that  vorticity  generation  occurs  in  the  flow  field  as  explained  by  the  fact 
that  the  equation  of  vorticity  transport  has  a  source  term  proportional  to  37  /  dx,  the  physical 
origin  of  which  lies  in  the  baroclinicity  of  the  flowfield.  The  isotherms  in  the  channel  show 
a  strong  tendency  to  become  horizontal  and  the  temperature  increases  in  the  positive  y 
direction. 

At  t  =  0. 15,  there  is  a  net  downward  volume  flux,  Q,  in  the  channel.  Many  recirculating 
flow  structures  persist,  although  the  tendency  to  decay  by  viscous  action  is  evident.  The 
isotherms  in  the  center  of  the  slot  are  horizontal  and  increasing  in  the  positive  y  direction. 
It  is  observed  that  the  motion  pressure  contours  are  normal  to  the  plates  and  to  the  sidewalls 
of  the  enclosure,  i.e,  3P  /  dx  is  negligible  in  the  slot  and  near  the  sidewalls  of  the  enclosure. 
At  t  =  0.2,  the  net  volume  flux  is  upward,  and  attains  a  local  maximum,  but  the  other  trends 
described  for  the  previous  case  continue.  At  t  =  1 .0  the  solution  approaches  the  steady  state 
solution.  At  this  time,  viscous  action  results  in  the  decay  of  several  of  the  vortices  and  many 
of  the  recirculating  flow  structures  coalesce.  There  are  two  counter  rotating  vortices  near 
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Figure  3.23  v-velocity  profiles  along  a  horizontal  line  at  the  center  of  the 
domain,  i.e.,  along  y  =  1 .5.  The  dotted  line  shows  the  location  of  the  heated  plate. 
Only  the  left  half  of  the  domain  is  shown.  The  profiles  clearly  show  the 
existence  of  a  thin  wall  layer  adjacent  to  the  heated  channel  walls.  The  letters 
E  and  F  correspond  to  the  physical  locations  shown  in  Fig.  3.1. 
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the  top  wall  of  the  enclosure  and  these  extend  downward  into  the  slot  region.  It  is  also  noted 
that  the  motion  pressure  variation  in  the  x-direction  is  negligible. 

A  qualitative  comparison  of  the  current  results  can  be  made  against  the  results  of  Chang 
and  Lin  (1989).  They  have  presented  detailed  studies  of  the  flow  and  thermal  structures  for 
the  case  of  Pr  =  0.7,  Ra  =  104,  /R  =  5.  A  monotonic  variation  of  the  flow  characteristics  up 
to  steady  state,  especially  the  recirculating  cells  in  the  early  stages  of  the  transient,  was 
obtained,  as  opposed  to  the  current  results,  where,  an  oscillatory  volume  flux,  decaying  to 
steady  state  can  be  observed.  However,  in  the  heat  transfer  results,  notably,  the  Nusselt 
number,  they  recorded  an  undershoot  in  the  time  evolution  of  the  Nusselt  number.  In  the 
current  calculations,  the  Nusselt  number  decays  monotonically  to  the  steady  state  and  no 
undershoot  can  be  observed.  Also,  Chang  and  Lin  (1989)  did  not  obtain  recirculating  flow 
in  the  channel  at  steady  state,  probably  due  to  the  higher  aspect  ratio  channel  used  in  their 
calculations.  The  recirculating  flow  observed  at  steady  state,  in  the  present  case,  is 
qualitatively  consistent  with  the  results  obtained  by  Shyy  et  al  (1992)  for  their  case  of  Pr  = 
0.71,  Gr  =  105,  /R=  1.  Low  Ra  cases  have  been  presented  by  Ramanathan  and  Kumar 
(1991),  where,  no  backflow  was  found.  They  have  also  presented  a  case  for  Ra  =  105,  for 
which  they  have  not  presented  the  complete  details  of  the  flow  and  thermal  fields.  However, 
they  show  that  the  temperature  distribution  along  the  channel  walls  were  qualitatively 
consistent  with  the  similarity  solution  of  Sparrow  and  Gregg  (1956).  This  agreement  with 
a  theoretical  prediction  based  on  the  parabolic  equations  implies  the  absence  of  recirculating 
flow  in  the  channel,  in  their  case. 

3.4  Steady  Thermocapillary  Convection  in  a  Liquid  Bridge 

Apart  from  buoyancy,  the  other  important  physical  mechanism  found  in  many  crystal 
growth  processes  is  capillarity  resulting  from  a  free  surface.  The  formulation  aspects  of  the 
float  zone  free  surface  have  been  presented  in  chapter  2.  A  major  goal  of  the  present  thesis 
is  to  simulate  the  complete  crystal  growth  process  in  a  float  zone.  Such  a  task  will  be 
presented  in  chapter  5  after  the  phase  change  aspect  has  been  addressed.  In  the  following, 
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thermocapillary  convection  in  a  liquid  bridge  has  been  conducted.  Such  a  liquid  bridge 
illustrated  in  Fig.  3.24(a)  may  be  considered  as  a  model  for  the  study  of  convective  heat 
transfer  mechanisms  in  the  float  zone  crystal  growth  process.  The  most  influential  factors 
on  the  crystal  quality  and  process  economics  are  the  heat  and  solute  transport  processes  in 
the  melt,  which  control  the  morphology  and  movement  of  the  solidifying  front.  In  the 
following  study,  a  configuration  frequently  encountered  in  practical  float-zone  processes  is 
isolated  for  detailed  analysis  of  heat  transport  within  the  melt.  Of  particular  interest  is  the 
role  of  buoyancy-driven  and  thermocapillary  convection  processes  and  the  interaction  of 
these  processes  with  the  shape  of  the  free  surface  in  determining  the  heat  transfer  within  the 
melt.  The  Navier-Stokes  equations,  along  with  the  energy  equation  are  written  in  cylindrical 
coordinates  to  facilitate  the  treatment  of  the  axisymmetric  geometry  shown  in  Figure 
3.24(b).  The  numerical  scheme  involves  the  discretization  of  the  transformed  form  of  the 
governing  equations  based  on  a  control  volume  formulation  as  described  earlier  in  chapter 
2.  The  free  surface  shape  is  chosen  so  as  to  satisfy  the  static  force  balance  between  the 
hydrostatic  pressure  of  the  melt  and  surface  tension  forces  generated  due  to  the  presence  of 
the  free  surface  between  the  melt  and  the  ambient  gas.  The  details  of  the  integration  of  the 
Young-Laplace  equation  to  arrive  at  an  equilibrium  shape  of  the  free  surface  have  been 
considered  previously  in  chapter  2. 

3.4.1    Scaling.  Nondimensionalization  and  Governing  Equations 

An  appropriate  choice  of  scales  for  the  nondimensionalization  procedure  is  arrived  at 
by  selecting  the  base  radius  rb  as  the  characteristic  length  scale  and  the  thermal  diffusion 
velocity  scale,  defined  as  a/r^,  where  a  is  the  thermal  diffusivity.  The  selection  of  the 
temperature  scale  is  not  as  straightforward  because  the  boundary  condition  at  the  free  surface 
is  specified  as  a  heat  flux  from  which  an  estimate  of  the  maximum  temperature  in  the  domain 
has  to  be  derived.  The  procedure  is  explained  below  in  detail.  Since  the  heat  flux  is  specified 
on  the  free  surface,  a  temperature  scale  has  to  be  derived  in  order  to  nondimensionalize  the 
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energy  equation  and  to  define  the  Grashof  and  the  Marangoni  numbers.  It  would  seem 
natural  to  define  a  temperature  scale  based  on  the  heat  flux 


q 


AT  =  ^  rb  (3.12a) 

where  qmax  is  the  maximum  value  of  the  heat  flux  specified  on  the  boundary.  However,  it 
has  been  observed  that  such  a  definition  consistently  overestimates  the  temperature  scale, 
that  is,  the  maximum  temperatures  achieved  in  the  domain,  even  for  conduction  cases,  are 
of  the  order  of  five  times  less  than  the  value  suggested  by  Eq.  (3.12a).  In  addition,  it  is 
expected  that  the  maximum  temperatures  obtained  by  taking  convection  into  account  should 
be  even  less.  Therefore,  the  temperature  scale  obtained  above  will  be  scaled  down  by  a  factor 
of  ten  in  order  to  obtain  a  more  realistic  temperature  scale.  Thus,  the  temperature  scale  used 
to  calculate  the  Grashof  and  Marangoni  numbers  is: 

AT  =  0.1  x  SE2i  rb  (3.12b) 

Regarding  the  convection  mechanisms,  the  buoyancy  effect  within  the  melt  and  the 
thermocapillary  effect  on  the  free  surface  have  been  included  in  this  model.  The  strength 
of  buoyancy  induced  convection  is  indicated  by  the  Rayleigh  number,  defined  as 

Ra  =  Gr  •  Pr  (3.13a) 

where  Gr  is  the  Grashof  number, 

grjjAg 

Gr  =  — (3.2b) 

Invoking  the  Boussinesq  approximation,  this  becomes 

gPxATrJ 

Gr  =  (3.14) 

where  Pr  is  the  Prandtl  number,  as  defined  in  Eq.  (3.2c).  The  relative  magnitude  of  the 
thermocapillary  effect  is  described  by  the  Marangoni  number,  Eq.  (3.4).  The  governing 
equations  may  be  nondimensionalized  with  respect  to  the  reference  scales  developed  in  Eq. 
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(2.5)  as  follows: 

(i)  continuity: 

(ii)  radial  momentum: 


d(Ru)  d(Ruv) 


dR 


dY 


(iii)  axial  momentum: 
d(Ruv) 


dR 
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d(Rv2) 


dY 


dmy  +  dim  =  0 
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R  Ra  Pr  0 


(3.15a) 


-  Pr  ^  (3.15b) 


(3.15c) 


(iv)  energy: 


a(RG0)  a(Rv0) 


dR 


dY 


_d_ 

dR 


(3.15d) 


Here  u  and  v  are  the  non-dimensional  velocity  components  in  the  radial  and  axial  directions, 

respectively  and  the  reference  length  scale,  d  is  set  to  rt>.  The  Boussinesq  approximation  is 
employed;  thus,  density  variations  are  neglected  in  the  governing  equations,  except  for  the 
buoyancy  terms  in  the  momentum  equations.  Only  the  right  half  of  the  domain  is  considered 
in  the  calculation  procedure,  and  symmetry  boundary  conditions  are  applied  at  the 
centerline.  The  associated  boundary  conditions  are  listed  below  in  nondimensional  form: 

d\ 


symmetry: 


At  R  -  °-  0  =  1  =  f  - 0 


top  and  bottom  boundaries: 

At  Y  =  0  and  1  :  0  =  0,  u  =  v  =  0 

At  the  free  surface;  R  =  f(Y): 

(i)  heat  flux  through  the  free  surface  (dimensionless) : 

q  =  V0  •  n 


(3.16a) 


(3.16b) 


(3.16c) 


which  equals  the  normal  derivative  of  the  temperature  at  the  free  surface,  and  n  is  the  unit 
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normal  vector  to  the  free  surface. 

(ii)  Marangoni  effect:  V  (V  ■  T)  •  n  =  Ma  V0  •  T  (3.10) 

where  the  shear  stress  equals  the  surface  tension  gradient  and  gives  rise  to  the  normal 
derivative  of  the  tangential  velocity  at  the  free  surface.  V  is  the  velocity  vector,  n  is  the  unit 

normal  vector  to  the  free  surface,  and  t  is  the  unit  tangent  vector  to  the  free  surface. 

(iii)  kinematic  condition  at  the  free  surface: 

V  •  n  =  0  (3.16d) 

which  specifies  a  no-penetration  condition  at  the  free  surface. 

3.4.2  Heat  Transfer  Calculations 

A  series  of  calculations  was  carried  out  with  the  geometry  of  the  free  surface  as  shown 
in  Figure  3.24.  The  material  selected  for  simulation  is  NiAl  whose  thermophysical 
properties  of  this  material  are  listed  in  Table  3.1.  If  phase  change  effects  are  not  considered, 
and  the  power  input  from  the  heater  is  2.3  X  106  W/m2,  then  the  temperature  scale  can  be 
obtained  from  Eq.  (3. 12b)  as  AT  =  36.5  K.  For  the  length  scale  defined  by  rb  =  1 .25  X  10"2 
m  as  in  the  previous  calculations  for  the  meniscus  shape,  this  temperature  scale  implies  a 
Marangoni  number  of  103,  a  Grashof  number  of  520,  and  a  Prandtl  number  of  0.04.  These 
values  vary  with  the  composition  of  the  NiAl.  Nevertheless,  it  is  clear  that  for  the  length 
scales  under  consideration,  thermocapillary  convection  plays  a  dominant  role  in  the  heat 
transfer  process.  In  the  computations  to  follow,  a  range  of  nondimensional  parameters, 
0  <  Gr  <  106,  0  <  Ma  <  103,   and  Pr  =  0.04,  0.1 

have  been  considered  for  investigation. 

3.4.3  Numerical  Procedure 

The  governing  equations  (3.15)  are  solved  using  the  pressure-correction  formulation  in 
nonorthogonal  body-fitted  coordinates.  The  convection  and  diffusion  terms  in  both  the 
momentum  and  energy  equations  are  discretized  using  second  order  central  differences.  The 
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(a)  Schematic  of  the  float-zone  process.  The  melt  zone  is  sustained 
by  the  surface  tension  between  the  melt  and  the  encapsulant.  The 
heating  profile  is  schematically  illustrated  in  the  right. 
Axisymmetry  has  been  assumed. 


0  =  0 


(b)  Schematic  of  the  mathematical  model  and  geometry  employed 
in  the  calculation  procedure. 


Figure  3.24  Schematic  of  the  float-zone  process  and  the  mathematical 
model  employed. 
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grid  points  were  clustered  towards  the  free  surface  and  for  the  computations  involving 
Marangoni  convection,  they  also  were  clustered  towards  the  lower  boundary  at  Y  =  0.  For 
the  cases  involving  natural  convection  only,  two  grid  systems  involving  81x52  nodes  and 
161  x  103  nodes,  were  used  to  discretize  the  domain.  The  solutions  are  indistinguishable, 
indicating  grid  independence.  The  solutions  presented  are  based  on  the  grid  system 
involving  161  X  103  nodes  for  cases  not  involving  Marangoni  convection  and  401  X  201 
nodes  for  cases  involving  Marangoni  convection 


Table  3.1  Thermophysical  properties  of  NiAl  and  process  parameters 


Density 

6000  kg/mJ 

Thermal  conductivity 

80  W/m-K 

Heat  capacity 

586  J/kg-K 

Coefficient  of  thermal  expansion 

1.51  X  10~5  K-1 

Kinematic  viscosity 

9  X  10"7  m2/s 

Surface  tension  coefficient 

-  2.7  X  10"4  N/m-K 

Surface  tension 

1.35  N/m 

Latent  heat 

5.1  X  105J/kg 

Melting  point 

1900  K 

Gravitational  acceleration 

9.81  m/s2 

3.4.4  Heat  Conduction 

A  calculation  was  carried  out  with  the  heat  transfer  taking  place  through  conduction 
effects  only.  The  purpose  was  to  carry  out  an  a  posteriori  verification  of  the  temperature 
scales  used  to  define  the  Grashof  and  the  Marangoni  numbers  in  subsequent  calculations. 
It  may  be  noted  that  due  to  the  curved  geometry  of  the  free  surface,  the  location  of  the 
maximum  temperature  shifts  towards  the  convex  portion  of  the  free  surface,  as  is  evident 
from  the  region  labelled  B  in  Figure  3.25.  The  nondimensional  value  of  0max  obtained  here 
is  2.3  based  on  the  temperature  scale  defined  in  Eq.  (3.12b). 

3.4.5  Natural  Convection 

Figure  3.26  shows  the  results  of  incorporating  the  effects  of  buoyancy  induced 
convection  with  a  Grashof  number  of  106  and  a  Prandtl  number  of  0.1.  The  streamline 
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B 


Isotherms:  max:  2.3,  20  contours 


Figure  3.25.  Temperature  distribution  in  the  float  zone  neglecting  convective  heat 
transport. 


pattern  shows  a  single  recirculating  zone  in  each  half  of  the  domain  and  is  tangential  to  the 
free  surface.  It  may  be  recalled  that  since  surface  tension  effects  have  not  been  accounted 
for,  the  free  surface  does  not  exert  a  shear  stress  on  the  fluid.  The  fluid  adjacent  to  the  free 
surface  is  heated  and  rises  to  the  top  of  the  domain  due  to  buoyancy  effects,  causing  the 
isotherms  to  be  distorted  towards  the  top  of  the  domain  and  the  location  of  @max  to  shift 
correspondingly  upwards.  The  magnitude  of  0max  drops  to  1.14,  reflecting  the  validity  of 
the  temperature  scale  defined  by  Eq.  (3.12b). 
3.4.6    Interaction  of  Natural  and  Thermocapillary  Convection 

Gr  =  106.  Pr  =  0.1.  and  Ma  =  500.  Figure  3.27(a)  shows  the  streamfunction  and 
isotherms  by  imposing  a  Marangoni  number  of 500.  Since  the  surface  tension  decreases  with 
temperature,  the  tendency  of  the  surface  tension  gradient  induced  shear  stress  is  to  reduce 
the  convection  strength  in  the  vicinity  of  the  convex  part  of  the  free  surface  (labelled  B  in 
Fig.  3.27a)  and  increase  the  convection  strength  near  the  upper  part  of  the  domain  where  the 
free  surface  is  concave  (labelled  A  in  Fig.  3.27a).  The  overall  effect  is  to  marginally  shift 
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Gr  =  106,  Ma  =  0,  Streamfunction:  max:  6.31,  20  contours 


Gr  =  106,  Ma  =  0,  Isotherms:  max:  1.14,  20  contours 


Figure  3.26.  Effect  of  natural  convection  on  transport  characteristics  in  the  float 
zone. 


the  location  of  0max  in  the  downward  direction,  towards  the  convex  portion  of  the  free 
surface  as  compared  with  the  Ma  =  0  solution  shown  in  3.27(b). 

Gr  =  106.  Pr  =  0. 1.  and  Ma  =  1000.  Figure  3.27(b)  shows  the  effect  of  increasing  the 
Marangoni  number.  The  trends  established  in  the  previous  case  continue,  shifting  the 
location  of  0max  downwards  towards  the  convex  portion  of  the  free  surface.  Since  the 
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Streamfunction 

max:  5.99,  20  contours 


Isotherms 

max:  1.14,  20  contours 


Figure  3.27.  Effect  of  Grashof  number,  Marangoni  number  and  Prandtl  number 
on  transport  characteristics  in  the  float  zone. 

(a).  Gr  =  106,  Ma  =  500  and  Pr  =  0.1.  The  shear  stress  due  to  the  surface  tension 
gradient  marginally  decreases  the  overall  convection  strength  since  it  counteracts 
the  buoyancy  induced  convection. 
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Streamf unction 

max:  5.55,  20  contours 


Isotherms 

max:  1.18,  20  contours 


Figure  3.27  —  continued.  Effect  of  Grashof  number,  Marangoni  number  and 
Prandtl  number  on  transport  characteristics  in  the  float  zone. 

(b).  Gr  =  106,  Ma  =  1000  and  Pr  =  0.1.  The  trend  exhibited  in  case  (a)  continues. 
The  increased  shear  stresses  due  to  the  stronger  surface  tension  gradient  causes  a 
noticeable  decrease  in  the  convection  strength  and  an  increase  in  maximum 
temperature  on  the  free  surface. 
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convection  strength  in  the  vicinity  of  0max  is  further  reduced,  the  value  of  0max  also 
increases. 

Gr  =  2000.  Pr  =  0. 1.  and  Ma  =  1000.  Figure  3.27(c)  shows  the  effect  of  decreasing 
the  relative  strength  of  buoyancy  induced  convection  versus  Marangoni  convection.  In  this 
case,  the  convection  pattern  is  dominated  by  the  Marangoni  effect.  The  convection  pattern 
consists  of  counterclockwise  rotating  convection  roll  at  the  top  right  of  the  domain  adjacent 
to  the  free  surface  (labelled  A)  and  clockwise  rotating  convection  roll  at  the  bottom  right  of 
the  domain  corresponding  to  region  labelled  B.  The  convection  rolls  adjacent  to  the  convex 
portion  of  the  free  surface,  which  occur  at  the  lower  portion  of  the  boundary,  are  substantially 
stronger  because  of  the  stronger  temperature  gradients  along  the  convex  portion  of  the  free 
surface.  The  location  of  0max  is  close  to  that  of  the  pure  conduction  solution  as  compared 
to  the  high  Gr  cases.  The  effect  of  Marangoni  convection  is  to  shift  the  location  of  0max t0 
the  convex  portion  of  the  boundary. 

Gr  =  2000.  Pr  =  0.04.  and  Ma  =  1000.  Figure  3.27(d)  shows  the  effect  of  decreasing 
the  Prandtl  number.  The  effect  of  lower  Prandtl  number  is  to  increase  the  convective  effects 
in  the  flow  (Shyy  1994)  by  decreasing  the  magnitude  of  the  viscous  terms  in  the  momentum 
equations.  In  the  energy  equation,  this  effect  will  be  balanced  by  the  increased  thermal 
diffusivity.  The  convection  pattern  is  qualitatively  similar  to  the  case  shown  in  Fig.  3.27(c). 
It  can  be  observed  that  the  isotherms  are  distorted  near  the  convex  portion  of  the  free  surface, 
due  to  the  strong  fluid  convection,  whereas  in  the  upper  portion  of  the  domain,  the  pattern 
resembles  the  pure  conduction  case. 

From  the  above  cases,  some  useful  inferences  can  be  made  that  are  directly  applicable 
to  the  float-zone  growth  of  NiAl.  For  the  length  scales  of  current  experimental  interest,  it 
is  observed  that  the  pure  conduction  model  of  heat  transfer  is  inadequate.  Strong  convective 
heat  transfer  effects  can  be  observed,  which  are  dominated  by  the  Marangoni  effect.  The 
Marangoni  effect  causes  clockwise  recirculating  convection  rolls  in  the  bottom  right  and  top 
left  of  the  domain  and  counterclockwise  recirculating  rolls  at  the  top  right  and  bottom  left 
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Streamfunction 

min:  -1.72  (-),  20  contours 


Isotherms 

max:  1 .47,  20  contours 


Figure  3.27 — continued.  Effect  of  Grashof  number,  Marangoni  number  and  Prandtl 
number  on  transport  characteristics  in  the  float  zone. 

(c).  Gr  =  2000,  Ma  =  1000  and  Pr  =  0. 1 .  Decreasing  the  Grashof  number  results  in 
a  clockwise  rotating  convection  cell  in  the  convex  region  of  the  meniscus  labelled  B. 
The  location  of  the  maximum  temperature  shifts  downward.  The  convection  cells  in 
the  upper  region  labelled  A  is  substantially  weaker  than  the  cases  dominated  by 
buoyancy  induced  convection. 
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Streamfunction 

min:  -1.63,  20  contours 


Isotherms 

max:  1.57,  20  contours 


Figure  3.27 — continued.  Effect  of  Grashof  number,  Marangoni  number  and  Prandtl 
number  on  transport  characteristics  in  the  float  zone. 

(d)  Gr  =  2000,  Ma  =  1000  and  Pr  =  0.04.  With  decreasing  Prandtl  number,  the 
magnitude  of  the  maximum  temperature  increases,  but  the  convection  cells  do  not 
penetrate  as  deeply  into  the  melt  zone  and  the  overall  convection  strength  is  reduced. 
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of  the  domain.  For  the  high  Bond  numbers  that  prevail  under  1-g  conditions  on  earth,  the 
free  surface  has  a  convex  shape  in  the  lower  portion  of  the  melt  due  to  the  hydrostatic 
pressure  of  the  melt  column.  This  convex  shape  induces  asymmetry  in  the  shape,  spatial 
extent,  and  strength  of  the  convection  rolls  and  the  distribution  of  the  temperature.  For  the 
pure  conduction  case,  the  convex  shape  causes  the  location  of  the  maximum  temperature  to 
be  shifted  towards  the  convex  side  of  the  free  surface.  In  the  presence  of  Marangoni 
convection,  this  effect  is  augmented  by  the  strong  convective  effect  near  the  convex  portion 
of  the  domain,  due  to  the  stronger  temperature  gradients  along  the  convex  portion  of  the  free 
surface.  This  convection  in  turn  causes  the  location  of  the  maximum  temperature  to  move 
further  downwards.  Marangoni  convection  completely  overwhelms  the  buoyancy-driven 
convection  for  this  particular  configuration  and  material  properties,  augmenting  the 
buoyancy  effect  in  the  upper  portion  of  the  domain  and  counteracting  it  in  the  lower  portion 
of  the  domain.  Thus,  the  dependence  of  surface  tension  on  temperature  plays  a  dominant  role 
in  the  heat  transfer  characteristics  in  the  float  zone. 


It  may  be  expected  that  convection  in  the  melt  will  impact  the  meniscus  shape  through 
the  normal  stress  terms  which  will  be  balanced  by  the  free  surface  curvature.  In  addition, 
the  free  surface  will  no  longer  be  an  isotherm,  which  implies  that  the  surface  tension  will  vary 
from  point  to  point  depending  on  the  spatial  distribution  of  the  temperature  along  the  free 
surface.  Under  such  conditions,  Equation  (2. 15)  must  be  modified  as  follows: 


where  u  is  the  viscosity  of  the  melt.  This  can  be  nondimensionalized  in  a  manner  similar 
to  Eq.  (2. 15)  along  with  the  reference  velocity  scale  defined  as  a/rt,,  where  a  is  the  thermal 
diffusivity.  In  nondimensional  form,  this  becomes: 


3.5   Effect  of  Convection  on  Meniscus  Shape 


(1  +  f'2)V2      f  (i  +  f'2)l/2 


r  i 


(3.17) 


Bo  y  -  AP  +  Ca^r1  =  ( 1  - 


dY  AT 


f    1 


(3.18) 


dT  o0 


(1+f2)3/2      f  (1  +  f  2)l/2 
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where  Ca  is  the  Capillary  number  defined  as  (2  \i  a)/[y0  rb),  Yo  is  the  surface  tension  at  the 
reference  temperature,  and  the  italics  denote  the  nondimensionalized  variables.  The 
Capillary  number  and  the  Bond  number  determine  the  relative  influence  of  the  hydrostatic 
and  the  convective  effects,  respectively.  The  term  premultiplying  the  curvature  term 
describes  the  correction  due  to  the  variation  of  surface  tension  with  temperature.  In  many 
instances  as  will  be  illustrated  later  it  is  found  that  the  Capillary  number,  Ca,  as  well  as  the 
normal  gradient  of  the  normal  velocity  are  both  small  compared  to  the  hydrostatic  term,  from 
an  a  posteriori  estimation,  and  hence  meniscus  shapes  based  on  Eq.  (2. 15)  generally  suffice. 

The  deformation  of  the  free  surface  due  to  thermocapillary  convection  is  dominated  by 
the  Capillary  number.  For  example,  Zebib  et  al.  (1985)  found  that  whereas  the  shape  is 
qualitatively  sensitive  to  the  Prandtl  number,  the  magnitude  of  the  deformation  is 
determined,  to  the  leading  order,  by  the  Capillary  number.  More  insight  into  the  effects  of 
Marangoni  convection  on  transport  dynamics  can  be  found  in  Carpenter  and  Homsy  (1989) 
and  Mundrane  and  Zebib  (1994).  Equation  (3.18)  describes  the  role  of  the  normal  stress, 
arising  from  the  normal  derivative  of  the  normal  velocity  component  at  the  free  surface,  in 
changing  the  shape  of  the  meniscus.  Under  1-g  conditions  and  at  the  length  scales 
considered,  the  hydrostatic  terms  dictated  by  the  high  Bond  numbers  are  the  dominant  terms 
in  Eq.  (3.18)  and  the  normal  stress  terms  due  to  convection  are  small  compared  to  the 
hydrostatic  term.  As  an  example  consider  the  impact  of  convection  on  the  shape  of  the  liquid 
bridge  considered  earlier  in  this  chapter.  The  Bond  number  for  this  case  is  Bo  =  8.8  and  the 
Capillary  number  is  Ca  =  1.5  X  10-5.  The  impact  on  the  surface  curvature  varies  with  the 

C    dV  Id 

vertical  distance  and  is  described  by  1  +  \ — .  This  quantity  may  be  plotted  to  show 

Bo  y 

the  variation  of  its  magnitude  along  the  y-axis  as  shown  in  Figure  3.28.  This  shows  that 
except  at  the  lowermost  corners  of  the  domain,  the  normal  stress  term  is  vanishingly  small. 
Considering  that  the  meniscus  is  pinned  at  both  the  top  and  the  bottom  of  the  domain,  the 
correction  to  the  overall  shape  of  the  meniscus  is  expected  to  be  modest  except  at  the 
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Figure  3.28  Semilog  plot  of  Y  versus  1  +  ^a  dVn/dn 

Bo  y 

This  shows  the  locations  where  the  normal  stress  term  due  to  convection  becomes 

important  compared  to  the  hydrostatic  pressure. 

lowermost  corner  where  the  curvature  will  be  substantially  lower.  Since  the  effect  is 
extremely  localized  at  the  corner  locations  where  the  flow  turns  to  match  the  boundary,  the 
overall  effect  on  the  heat  transfer  and  convection  characteristics  is  expected  to  be  negligible 
for  this  case.  Accordingly,  corrections  to  the  meniscus  shape  appear  unnecessary  in  this 
particular  study.  We  also  observe  that  the  percentage  variation  of  the  surface  tension  with 
temperature  is  less  than  1%  and  hence  cannot  significantly  affect  the  force  balance  at  the 
interface. 

Similar  results  have  been  reported  by  Tao  et  al.  (1995),  where  a  vorticity-stream 
function  formulation  was  employed.  The  fluid  chosen  in  that  study  is  silicone  oil,  which  has 
a  very  different  Prandtl  number,  of  the  order  of  102  to  104.  Their  solutions  can  help 
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understand  Marangoni  convection  in  a  different  parameter  range.  In  a  recent  effort,  Hou  et 
al.  (1994)  have  employed  a  boundary  integral  approach  to  simulate  surface  tension  driven 
flows,  with  good  results  obtained  for  a  wide  range  of  phenomena.  Kanouff  and  Greif  (1994) 
have  analyzed  oscillations  in  thermocapillary  convection  in  a  square  cavity. 

3.6  Concluding  Remarks 

In  this  chapter,  different  convection  schemes  have  been  investigated  and  applied  in  the 
context  of  the  SIMPLE  algorithm  to  compute  fluid  convection  and  heat  transfer 
characteristics.  Two  test  cases  have  been  considered  to  illustrate  the  computational  and 
physical  aspects  of  buoyancy  driven  and  thermocapillary  convection.  The  problem  of 
electronics  cooling  in  a  cabinet  has  been  modeled  by  considering  transient  natural 
convection  inside  an  enclosed  vertical  channel.  Accurate  simulations  have  been  conducted 
to  capture  the  transient  effects  which  are  important  during  startup  operations.  The  second 
problem  considered  is  Marangoni  convection  in  a  liquid  bridge  serving  as  a  simplified  model 
of  the  float  zone  crystal  growth  process.  These  simplified  computations  serve  as  the 
foundation  for  more  elaborate  and  complex  float  zone  simulations  later  in  this  study. 

This  chapter  also  addresses  the  effect  of  melt  convection  on  the  free  surface  shape  in  the 
context  of  small  Capillary  numbers  encountered  in  the  float  zone  growth  process.  It  is  found 
that  in  many  practical  configurations  under  1-g  conditions,  the  combination  of  high  Bond 
number  and  small  Capillary  numbers  results  in  negligible  impact  of  melt  dynamics  on  the 
free  surface  shape.  These  calculations  form  the  building  block  for  the  development  of  a  more 
complete  algorithm  for  float  zone  simulation. 


CHAPTER  4 

ENTHALPY  BASED  PHASE  CHANGE  MODELING 


4.1  Background 

Apart  from  establishing  the  capability  of  solving  the  fluid  flow  and  convective  heat 
transfer  problem,  the  computation  of  phase  change  is  the  other  major  challenge.  In  this 
chapter,  a  numerical  method  is  investigated  and  enhanced  to  compute  phase  change 
phenomena  on  a  fixed  grid  in  the  context  of  pressure  based  algorithms.  In  this  method,  also 
known  as  the  enthalpy  formulation,  the  energy  equation  is  initially  cast  in  terms  of  the  total 
enthalpy.  By  splitting  the  total  enthalpy  into  the  sensible  enthalpy  and  the  latent  heat,  a 
unified  formulation  can  describe  energy  transport  in  both  phases.  The  relationship  between 
the  sensible  heat  and  the  latent  heat  is  modeled  in  a  simplified  manner  without  recourse  to 
either  morphological  or  microscopic  scales.  In  this  manner,  phase  boundaries  are  captured 
naturally  rather  than  by  separate  explicit  tracking  procedures.  The  spatial  averaging  in  the 
enthalpy  method  implies  that  a  phase  fraction  is  defined  based  on  the  volume  fraction 
occupied  by  a  given  phase,  say,  liquid.  Such  a  two-phase  approach  is  naturally  suitable  for 
alloys  where  phase  change  generally  occurs  over  a  temperature  range,  according  to  the  local 
mixture  composition,  forming  a  mushy  zone  in  which  both  phases  can  coexist.  A  detailed 
discussion  can  be  found  in  Crank  (1984).  For  pure  materials,  which  undergo  phase  change 
at  a  given  temperature,  the  mushy  zone  is  a  modeling  artifact.  This  practice  is  akin  to  the 
shock  capturing  technique  for  compressible  flows  where  a  finite  thickness,  larger  than 
physical  reality,  of  the  order  of  one  or  two  grid  spacings  is  needed  to  define  the  shock. 
However,  for  both  shock  and  phase  change  calculations,  this  artificial  smearing  is  generally 
quite  acceptable  and  does  not  necessarily  impact  the  global  accuracy  of  the  flow  solution. 
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A  major  requirement  for  both  problems  is  to  devise  numerical  schemes  that  are  accurate  and 
robust  in  the  bulk  of  the  flow  domain  and  can  resolve  the  internal  boundaries  (shock, 
interface  etc.)  within  one  or  two  computational  cells.  As  will  be  demonstrated  later,  this 
requirement  can  be  quite  satisfactorily  met  by  the  enthalpy  formulation. 

4.2  Heat  Conduction 
In  the  following,  the  enthalpy  formulation  and  the  associated  numerical  algorithm  is 
developed  via  application  to  a  heat  conduction  problem.  The  governing  equation  can  be  cast 
in  the  form 

=  V  •  (kVT)  (4.1) 

where  H  is  the  total  enthalpy.  The  equation  can  be  discretized  according  to  the  conventional 

control  volume  procedure,  wherein  the  total  enthalpy  H  and  the  temperature  T  can  be 

interpreted  as  averaged  values  within  a  control  volume.  Then  H  can  be  written  as 

H   =  Cp  T  +  f  •  L  (4.2) 

where  f  is  the  phase  fraction  of  phase  1  and  can  be  interpreted  as  the  fractional  volume  of 

a  computational  cell  occupied  by  that  phase.  The  phase  fraction  of  phase  1  is  zero  in  the 

region  occupied  by  phase  2  and  unity  in  the  region  occupied  by  phase  1 .  The  phase  fraction 

lies  between  zero  and  unity  when  the  control  volume  is  undergoing  phase  change.  L  is  the 

latent  heat.  Substituting  equation  (4.2)  into  the  governing  equation  (4.1)  we  get: 

d(pCpT)  d(pf) 

V       =  V  •  (kVT)   -  L  •  (4.3) 

where  the  latent  heat  now  appears  as  a  source  term.  This  formulation  is  a  single-region 
formulation,  wherein  one  set  of  governing  equations  can  describe  both  phases  and  hence  can 
more  readily  handle  complex  interface  shapes,  branched  and  multiple-valued  interfaces.  To 
close  the  physical  model,  a  unique  relationship  between  the  fluid  fraction  f  and  the 
temperature  T  has  to  be  formulated.  For  a  pure  substance  undergoing  isothermal  phase 
change,  the  total  enthalpy,  H,  is  a  discontinuous  function  of  the  temperature  (Fig.  4.1a) 
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H  = 


(4.4) 


[Cp  T  ,   T  <  Tm 
|Cp  T  +  L  ,   T  >  Tm 

However,  from  a  computational  viewpoint  discontinuities  are  difficult  to  track,  and  often  it 
is  necessary  to  smear  the  phase  change  over  a  small  temperature  range  to  attain  numerical 
stability,  as  schematically  illustrated  in  Fig.  4.1(b).  Thus, 

fCp  T  ,  T  <  Ts 


(T  -  Ts  )  ,   Ts  <  T  <  T, 


(4.5) 


H  =  <!Cp  Ts  + 

Cp  T  +  L  ,   T  >  T, 

where  Ts  =  Tm  -  e  ,  T\  =  Tm  +  e,  and  8  is  a  phase  change  interval.  The  discontinuity  is 
thus  replaced  by  a  small  interval  over  which  phase  change  occurs.  For  a  pure  substance 
undergoing  isothermal  phase  change,  this  is  a  purely  numerical  artifact  and  needs  to  be  as 
low  as  possible  in  order  to  accurately  model  the  physical  system.  The  above  relationship  is 
substituted  into  Eq.  (4.5)  to  obtain: 

T   -  T, 


f  = 


2  e 


(4.6) 


which  is  used  to  iteratively  update  the  phase  fraction  from  the  computed  temperature  field. 
This  technique  has  been  used  successfully  by  many  researchers,  including  Bennon  and 
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T  (b) 

T 

Figure  4. 1    (a)  Schematic  of  the  enthalpy  —  temperature  relationship  for  a  pure 
substance;  (b)  Illustration  of  the  numerical  treatment  of  the  discontinuity  which  is 
accomplished  by  smearing  the  phase  change  interval  over  a  finite  temperature  range. 
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Incropera  (1987),  Lacroix  (1989),  Lacroix  and  Voller  (1990),  Prakash  (1990),  Prakash  and 
Voller  (1989),  Shyy  and  Chen  (1991),  Voller  and  Cross  (1981),  Voller  and  Prakash  (1987) 
and  henceforth  will  be  referred  to  as  the  T-based  update  method.  An  alternative  formulation 
is  to  express  the  phase  fraction  as  a  function  of  the  total  enthalpy  rather  than  the  temperature, 
whereby  a  continuous  relationship  can  be  obtained  (Atthey  1975).  An  implicit  procedure 
is  examined  here.  First,  we  invert  the  H  =  H(T)  relationship  (4.5)  (see  Fig.  4.2)  so  that 

H 


,   H  <  Hs 

Ts  + 
H  - 


2e 


L  +  2Cpe 


(H  —  Hs)  ,   Hs  <  H  <  H, 


(4.7) 


H  >  H 


where  Hs  and  Hi  are  the  enthalpy  values  corresponding  to  the  Ts  and  the  Tj  temperatures, 
respectively,  Hs  =  Cp  Ts  and  Hi  =  Cp  Ti  +  L,  for  constant  Cp_  Substituting  the  above 
relationships  into  Eq.  (4.7),  we  obtain 

for  the  iterative  update  of  the  fluid  fraction.  The  motivation  for  the  above  formulation  stems 
from  the  fact  that  the  temperature  is  always  a  continuous  function  of  the  total  enthalpy.  In 


T 

1 

/ 

H 

Figure  4.2.    Sketch  showing  the  inverted  enthalpy  —  temperature 
relationship  for  a  pure  material. 
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particular,  it  is  noted  that  e  =  0  is  allowed  in  this  formulation,  thereby  accurately  modeling 
phase  change  of  a  pure  substance.  Another  motivation  is  the  fact  that  in  the  phase  change 
zone,  the  phase  fraction  is  a  rapidly  varying  function  of  the  temperature,  causing  it  to  be 
sensitive  to  small  errors  in  temperature.  In  particular,  this  can  result  in  an  unphysically  thick 
phase  change  zone  (mushy  region)  due  to  round-off  errors  in  the  temperature  calculation. 
Contrast  this  with  the  fact  that  if  the  latent  heat  is  large  relative  to  the  sensible  enthalpy,  the 
phase  fraction  is  a  slowly  varying  function  of  the  total  enthalpy  resulting  in  a  more  stable 
update  procedure.  This  update  procedure  will  henceforth  be  referred  to  as  the  H-based 
update  method. 

4.3  Numerical  Implementation 

4.3. 1    The  T-Based  Method 

A  straightforward  implementation  of  the  T-based  method,  using  the  notation  due  to 
Patankar  (1980),  would  be 

aPTp  =  £anbTnb  +  bk  —  q  •  L  [fir1  "  r1]  (4-9) 

followed  by  the  iterative  update 


T£  -Ts 


TP  <  Ts 


Ts  <  T£  <  T,  (4.10) 


2  e  '  ls-lp-1 
1  ,     T£  >  T, 


where  the  superscript  k  stands  for  the  current  iteration  value  at  a  given  time  level  n  and 
superscript  n  -  1  stands  for  the  previous  time  level.  For  Stefan  numbers  less  than  unity  (large 
values  of  latent  heat  relative  to  the  sensible  heat  CpAT),  this  form  of  the  iterative  update  is 
prone  to  numerical  instability  in  the  form  of  oscillations,  although  for  moderate  Stefan 
numbers,  underrelaxation  of  the  source  term  or  the  phase  fraction  update  facilitates  the 
solution  process.  Similar  findings  have  been  reported  by  Prakash  et  al.  (1987),  Prakash  and 
Voller  (1989),  Voller  and  Prakash  (1987).  The  chief  cause  of  this  oscillatory,  nonconvergent 
behavior  is  the  extreme  sensitivity  of  the  phase  fraction  to  small  changes  in  temperature, 
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leading  to  a  large  negative  feedback  effect  through  the  source  term.  This  tendency  gets  worse 
as  the  Stefan  number  is  reduced  and  the  time  step  is  refined.  A  simple  hypothetical  example 
will  serve  to  illustrate  the  deficiencies  of  the  above  procedure  and  also  will  provide  a 
criterion  for  devising  a  stable  update  procedure.  Consider  some  typical  nondimensional 
parameters:  L  =  20, Tm  =  0.5,  Cp=l.  Choose  e  =  0. 1 ,  then Ts  =  0.4  and  Ti  =  0.6.  Suppose 
that  in  the  course  of  the  calculation  we  obtain  Tk  =  0.39 ,  ft  =  0  and  Tk+1  =  0.42;  now  calculate 
the  total  enthalpy,  H  =  TK+1  +  ft  AH  =  0.42;  carrying  out  the  update  procedure,  we  get  ft+i 
=  0.1.  The  total  enthalpy  H  jumps  to  2.42,  an  increase  of  about  500%  due  to  the  update 
procedure  alone.  This  large  jump  in  total  enthalpy  due  to  the  update  procedure  is 
nonphysical;  the  temperature  and  fluid  fraction  fields  are  no  longer  consistent  with  the 
discrete  form  of  the  conservation  law.  Thus  this  update  procedure  will  not  be  able  to  yield 
mutually  consistent  T  and  f  fields  and  hence  is  not  guaranteed  to  converge  (Shyy  and  Rao 
1994a).  A  more  stable  update  procedure  can  now  be  devised  by  substituting  Eq.  (4. 10)  into 
the  discrete  form  of  the  conservation  law  and  moving  the  resulting  coefficients  of  Tp  into 
the  ap  term  to  give 


Tp"1  <  Ts 


apTP  =  XanbTnb  +  b*-o.L^  [f-" 
h  +  *  •  ^  4]  ^  "  I*»A  +  *  •  L^*  [J  +  J"']  .  (4.11) 


apTP=XanbTkb  +  b*-o.L^  [f^  -  1 


T,  <  Tl~l  <  T 


Tk_1  >  T, 


At 

where  the  superscript  k  indicates  the  current  iterative  value.  It  can  be  observed  that  with  this 
form  of  the  update  procedure,  the  temperature  and  the  fluid  fraction  field  are  computed 
simultaneously  and  are  therefore  mutually  consistent  at  every  step  of  the  iterative  process. 
This  procedure  converges  quickly  while  still  retaining  accuracy,  as  will  be  demonstrated 
later. 
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4.3.2    The  H-Based  Method 

The  implementation  of  the  H-based  method  is  relatively  straightforward.  Equation  (4.9) 
is  applied  at  every  discrete  cell,  followed  by  the  update  procedure,  which  can  be  written  as 
follows: 

(i)  Compute  the  total  enthalpy 

H£   =  CpT£  +   fjT1  •  L  (4.12) 

(ii)  Carry  out  the  update  procedure: 

0  ,     H*  <  Hs 

H£  -  Hs 


fk  =  J 


h,-h7  '    Hs  -  hp  -  Hi  (413) 

1  ,     H£  >  H, 


4.4  Results  and  Discussion 

4.4. 1    Accuracy  Assessment 

A  one-dimensional  conduction  driven  phase  change  problem  with  an  exact  solution  was 
chosen  as  the  test  bed  for  the  assessment  of  accuracy  (Chandra  1990).  A  schematic  is  shown 
in  Fig.  4.3.  Consider  an  idealized  phase  change  material,  with  constant  thermophysical 
properties,  including  density,  and  with  phase  change  occurring  at  a  constant,  specified 
temperature.  This  implies  the  absence  of  convection,  the  phase  change  being  driven  entirely 
by  the  conduction  process  and  release  of  latent  heat  at  the  phase  change  interface.  This  class 
of  problems  is  also  known  as  the  Stefan  problem.  The  material  is  initially  at  the 
phase-change  temperature,  Tm  =  0.  At  time,  t  =  0,  the  temperature  at  x  =  0  is  instantaneously 
raised  to  Tj,  =  1.  The  interface  is  defined  by  the  T  =  Tm  =  0  isotherm.  For  the  pure 
conduction  problem  in  each  phase  we  have 

=  ttj  V2Tj         i  =  liquid,  solid  (4.14a) 


along  with  the  jump  condition: 


,   dTs        ,  3Tj 
s~dn     ~     l~dn     =   Q  L  Vn  (4.14b) 
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Figure  4.3.  Schematic  of  the  phase  change  problem  showing  the  domain  and 
boundary  conditions.  For  the  conduction  driven  problem,  the  interface  is  planar. 
For  the  problem  with  gravity  driven  convection  included,  the  gravity  vector  is  in 
the  negative  y  direction. 


where  vn  is  the  interface  velocity.  For  this  problem,  a  closed-form  solution  exists  for  the 
time-dependent  location  of  the  phase  change  front  given  by  x  =  2X  /at,  where  \  is  the 


root  of 


A  exp(A2)  erf  (A)  =  (4.15) 

v   '  471 

and  St  is  the  Stefan  number;  defined  as  St  =  Cp  (Th  -  Tm)/L.  Note  that  the  heat  flux  to 
the  interface  takes  place  only  through  the  liquid  and  the  solid  remains  at  a  uniform 
temperature. 

Calculations  were  carried  out  on  uniform  grids  consisting  of  21  and  41  nodes  with  a 
time  step  At  =  0.01  and  on  a  grid  with  81  points  with  time  steps  At  =  0.1,  0.01,  and  0.001 
for  both  the  H-based  and  the  T-based  methods  in  order  to  compare  the  two  methods  on  the 
basis  of  accuracy.  Comparisons  were  made  with  the  exact  solution  to  assess  the  accuracy 
of  the  two  methods  on  the  basis  of  the  interface  location.  The  interface  location  was 
obtained  by  interpolating  for  the  f  =  0.5  contour.  For  the  T-based  method,  e  =  10~5  was 
adopted.  Second-order  central  differences  were  used  for  all  spatial  derivative  terms  and 
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first-order  Euler  backward  differencing  was  employed  for  marching  in  time.  Figures  4.4(a) 
and  4.4(b)  show  that  the  H-based  and  T-based  methods  yield  identical  solutions  even  for 
coarse  grids.  An  error  distribution  was  formed  by  computing  the  relative  error  at  every  time 
instant.  Figure  4.4(c)  shows  that  the  error  distribution  oscillates  about  a  zero  mean.  This 
can  be  correlated  with  the  fact  that  a  significant  part  of  the  latent  heat  content  in  a 
computational  cell  is  released  when  it  starts  to  undergo  phase  change.  Thus  the  interface  first 
slows  down  when  a  computational  cell  just  starts  to  change  state  and  then  subsequently 
speeds  up.  The  computed  solution  coincides  with  the  exact  solution  when  the  interface  is 
at  the  boundary  between  two  control  volumes.  The  amplitude  of  the  error  distribution  drops 
with  increasing  spatial  resolution,  and  its  frequency  doubles  when  the  grid  spacing  is  halved. 
A  measure  of  the  cumulative  error  is  obtained  from  the  norm  of  the  error  distribution,  which 
is  the  rms  value  of  the  error  distribution.  A  log-log  plot  of  the  cumulative  error  versus  Ax 
(Fig.  4.4(d))  shows  that  the  cumulative  error  varies  as  (Ax)1-5.  A  calculation  conducted  with 
a  time  step  of  0.001  yielded  a  solution  that  could  not  be  distinguished  from  the  solution  for 
the  time  step  of  0.01.  This  indicates  that  the  solution  is  time-accurate. 
4.4.2    Performance  Assessment 

The  problem  described  above  was  generalized  in  order  to  make  an  extensive  comparison 
of  the  two  update  procedures  formulated  earlier.  In  general,  the  solid  and  liquid  phases  differ 
in  their  thermal  diffusivities.  Thus  the  above  problem  was  reformulated  such  that  a  steady 
state  solution  would  exist  with  the  interface  positioned  within  the  computational  domain  so 
that  the  effects  of  heat  conduction  through  the  solid  would  also  play  an  equally  important 
role.  So  the  phase-change  temperature  was  set  to  Tm  =  0.5  with  all  other  parameters  being 
the  same  as  the  previous  study. 

Figures  4.5,  4.6,  and  4.7  show  the  solutions  and  the  relative  performance  obtained  for 
various  combinations  of  e,  At,  and  Ax.  Figures  4.5  and  4.6  show  the  solutions  and 
performance  comparisons  on  a  coarse  grid  of  Ax  =  0.04  for  time  steps  At  of  0. 1  and  0.01  and 
with  e  of  10-2,  10-3,  10-4,  and  10-6.  As  e  decreases,  the  step-like  character  of  the  interface 
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movement  becomes  obvious.  Corresponding  to  these  step-like  movements,  there  are  sharp 
spikes  in  the  iteration  count  required  for  convergence  at  the  given  time  step.  It  is  evident 
that  the  T-Based  method  requires  far  fewer  iterations  when  £  and  At  are  relatively  large.  It 
is  also  noted  from  Fig.  4.5(d)  that  the  T-Based  method  performs  better  than  the  H-Based 
method  as  the  interface  slows  down  and  approaches  the  steady  state.  For  a  combination  of 
low  e  and  smaller  At,  the  T-Based  method  does  not  converge  and  solutions  can  be  obtained 
only  with  the  H-Based  method  as  is  evident  in  Fig.  4.5(e)  and  (f)  and  also  in  Fig.  4.6(e)  and 
(f),  where  e  =  10-6  and  a  time  step  At  =  0.01  have  been  employed.  Figure  4.7  shows  the  effect 
of  grid  refinement.  The  spatial  resolution  was  doubled  to  Ax  =  0.02.  In  Fig.  4.7(a),  for  e 
=  10~2,  the  step-like  behavior  almost  vanishes  but  reappears  when  £  is  decreased  to  10-3  in 
(c).  The  corresponding  performance  curves  in  (b)  and  (d)  show  the  same  trends,  that  is,  the 
T-Based  method  performs  better  under  more  benign  conditions,  but  when  E  is  reduced  the 
relative  advantage  reduces  and  the  H-Based  method  becomes  more  efficient.  The  H-Based 
method  is  not  very  sensitive  to  the  value  of  E  and  can  work  for  E  =  0.  Also,  the  T-Based 
method  is  more  efficient  when  the  interface  motion  approaches  the  steady  state.  In  general, 
it  appears  that  the  T-Based  method  is  relatively  efficient  compared  to  the  H-Based  method 

L  Ax  Av    Tn        T'n —  1 

when  the  magnitude  of  the  nonlinear  term   —   2  £          *S  su^"lc^ent^y  ^ow-  ^n 

particular,  the  performance  degrades  rapidly  with  decreasing  £  until  a  point  is  reached  where 
a  solution  is  no  longer  feasible  with  the  T-Based  method. 

4.4.3  Summary 

The  enthalpy  formulation  has  been  applied  to  a  one-dimensional  heat  conduction  driven 
phase-change  problem  to  examine  alternative  procedures  for  the  phase  fraction  update.  The 
range  of  applicability  and  relative  performance  of  the  T-based  and  H-based  updates  have 
been  examined  in  the  context  of  solidification  problems.  It  is  shown  that  the  T-based  update 
can  yield  better  performance  provided  the  time  step  is  not  too  small  and  the  grid  has  sufficient 
resolution.  Recently,  Ding  and  Anghaie  (1994)  have  found  that  the  H-based  update  yields 
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Figure  4.4.  Solution  features  for  the  pure  conduction  problem  with  St  =  0.042:  (a) 
Coarse  grid  solution  (21  grid  points);  (b)  Fine  grid  solution  (41  grid  points).  In  each 
case  the  solid  line  designates  the  solution  obtained  using  the  H-based  method,  the 
dashed  line  using  the  T-based  method  and  the  dotted  line  is  the  exact  solution.  The 
H-based  and  the  T-based  method  yield  identical  results,  (c)  Error  at  every  time  instant 
for  various  grid  sizes.  The  solid  line  is  the  error  with  Ax  =  1/80,  the  dashed  line  is  the 
error  with  Ax=  1/40,  and  the  dotted  line  is  the  error  with  Ax  =  1/20.  The  error  oscillates 
about  a  zero  mean,  its  amplitude  decreases  and  the  frequency  doubles  when  the  grid 
spacing  is  halved,  (d)  Cumulative  error  versus  grid  spacing  for  At  =  0.01.  The  solid 
line  is  the  numerical  result,  and  the  dotted  line  has  an  approximate  slope  of  1.5. 
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Ax  =  0.04 
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T-Based 


(a)   e=  10  ,  At  =  0.1 


(b)  e=  10  ,  At  =  0.1 


(c)  e=10  ,  At  =  0.1 


(d)  e=  10  ,  At  =  0.1 


(e)  e=  10  ,  At  =  0.01 


(f)  e=  10  ,  At  =  0.01 


Figure  4.5.  Comparison  of  the  H-Based  and  the  T-Based  update  methods  for  the 
conduction  problem  using  a  grid  with  26  uniformly  distributed  points,  i.e.,  Ax  = 
0.04.  In  all  cases  the  solid  line  represents  the  H-Based  method  and  the  dash-dotted 
line,  the  T-Based  method,  (a)  and  (b)  show  the  solutions  obtained  with  e  of  IQr2 
and  10-3,  respectively,  (c)  and  (d)  show  the  relative  performance  of  the  two 
schemes  with  e  of  10~2  and  10~3,  respectively.  In  (e)  and  (f),  where  the  time  step 
has  been  reduced  from  At  =  0. 1  to  At  =  0.01,  the  T-Based  method  could  not  reach 
a  residual  of  10^*  and  has  been  excluded  from  consideration. 
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Figure  4.6.  In  (a)  and  (b)  the  value  of  e  has  been  further  reduced  to  10-4.  It  is 
seen  that  the  T-Based  method  does  not  yield  solutions  when  the  time  step  is 
reduced  from  At  =  1 .0  to  At  =  0. 1 .  When  the  T-Based  method  works,  it  gives  better 
performance  than  the  H-Based  method,  as  shown  in  (c)  and  (d).  However,  from 
(a)  and  (b)  it  is  obvious  that  the  solutions  are  not  time  step  independent,  (e)  shows 
the  results  of  the  H-Based  method  for  e  =  1 0~6  and  At  =  0.01.  The  T-Based  method 
does  not  converge  to  a  residual  level  of  10"4  for  this  combination  of  lower  e  and 
lower  At. 
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Figure  4.7.  This  shows  the  effect  of  grid  refinement.  The  calculations  were  carried 
out  on  a  grid  with  5 1  uniformly  distributed  points,  (a)  and  (c)  show  that  the  step-like 
behavior  of  the  solutions  decreases  with  increasing  spatial  resolution.  The  spikes  in 
the  performance  curves  in  (b)  and  (d)  correspond  to  the  step-like  movement  of  the 
interface  in  (a)  and  (c),  respectively.  The  T-Based  method  performs  better  in  both 
cases,  but  its  comparative  advantage  reduces  as  e  is  reduced  from  10~2  to  10~3.  Also, 
it  may  be  observed  that  the  H-Based  method  performs  better  when  the  temporal 
gradients  are  high,  i.e.,  when  the  interface  moves  faster  due  to  higher  thermal  fluxes. 
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better  performance  for  the  evaporation  problem  where  the  boiling  point  varies  with  the 
pressure  according  to  the  Clausius-Clayperon  equation  (Shyy  1994).  Thus,  the  choice  of  the 
update  procedure  depends  on  the  computational  parameters,  such  as  time  step  and  grid 
resolution,  and  on  physical  parameters  such  as  the  phase  change  interval  and  whether  the 
interface  temperature  varies  with  the  pressure,  as  in  the  evaporation/condensation  problem. 

4.5  Convective  Effects 

So  far,  only  conduction-dominated  effects  have  been  taken  into  account,  but,  as  will  be 
demonstrated  in  the  following,  convection  can  substantially  modify  the  flow  and  heat 
transfer  characteristics.  Consequently,  the  interface  shape  and  translation  velocity  can  be 
substantially  altered  by  fluid  convection.  Consideration  of  fluid  convection  implies  that  the 
fluid  fraction  field  be  coupled  with  the  flow  solver  along  with  the  energy  equation.  In  the 
energy  equation,  apart  from  the  fluid  convection  terms,  additional  source  terms  are 
incorporated.  In  the  momentum  equations,  the  velocities  in  the  solid  need  to  be  set  to  the 
casting  rate  of  the  solid  phase.  One  way  to  accomplish  this  is  by  considering  the  mushy  zone 
to  be  a  porous  medium  and  incorporating  D'  arcy  law  terms  in  the  momentum  equations.  The 
D'arcy  law  terms  vary  the  porosity  as  a  function  of  the  fluid  fraction  from  unity  in  the  liquid 
to  zero  in  the  solid. 
4.5.1    Governing  Equations 

The  governing  equations  in  Cartesian  coordinates  and  in  dimensional  form  for 
two-dimensional,  incompressible  flow  can  be  written  as: 

.     .  d(pu)  d(pv) 

continuity:  +  =  0  (4.16a) 


x-momentum: 
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y-momentum: 
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(4.16c) 
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(4.16d) 


In  the  above,  q  is  the  density,  u  is  the  velocity  component  in  the  x-direction,  v,  the  velocity 
component  in  the  y-direction,  p  is  the  pressure,  u.  is  the  viscosity,  k  is  the  conductivity,  g  is 
acceleration  due  to  gravity,  L  is  the  latent  heat  of  fusion,  Cp  is  the  specific  heat,  f  is  the  liquid 
fraction,  ucast  and  vcay,  denote  the  prescribed  velocity  at  which  the  solid  phase  is  moving, 
say,  for  example,  in  a  crystal  growth  or  a  casting  process.   Here,  the  gravity  vector 

A 

g  =  —  g  j  is  assumed  to  be  in  the  negative  y  direction.  For  inclined  configurations,  the 

appropriate  components  should  appear  in  both  the  momentum  equations. 
4.5.2    Source  Terms  in  the  Momentum  Equations 

When  computing  the  solidification  process  on  a  fixed  grid,  special  treatment  is  needed 
to  enforce  the  solid  velocity  to  be  either  zero  or  equal  to  the  pulling  rate  of  the  solid  phase 
in  a  typical  crystal  growth  or  casting  operation.  The  approach  taken  here  is  the  inclusion  of 
D'arcy-type  source  terms  in  the  momentum  equations.  Here  the  phase  change  material  is 
considered  to  be  a  porous  medium,  with  the  porosity  changing  from  zero  to  unity  as  the 
material  melts.  For  pure  materials,  the  porosity  changes  abruptly  as  the  phase  change  occurs, 
but  to  ensure  numerical  stability  a  continuous  variation  should  be  imposed.  Following  Shyy 
and  Chen  (1991a),  the  term  A  in  Eqs.  4. 16(b,c)  is  a  function  of  the  liquid  fraction  as  follows: 
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(4.17) 


where  the  constants  C  and  q  are  adjusted  so  that  A  is  at  least  seven  orders  of  magnitude  higher 
than  all  the  other  terms  in  the  momentum  equations  in  the  solid  region.  For  a  pure  material, 
which  melts  isothermally,  the  concept  of  a  porous  mushy  zone  is  a  numerical  artifact,  but 
for  alloys,  which  melt  over  a  temperature  range,  the  porous  mushy  zone  is  a  physical  reality, 
and  convection  within  the  mushy  zone  definitely  influences  the  structure  of  the  mushy  zone 
(Voller  and  Prakash  1987).  However,  it  is  not  straightforward  to  obtain  a  relationship 
between  the  morphology  of  the  mushy  zone  and  the  D'arcy  source  term,  (4.17),  and  it  is 
certainly  a  matter  for  further  investigation. 
4.5.3    Sources  of  Convection 

Depending  on  the  geometry  and  the  configuration,  there  can  be  several  driving  forces 
for  convection  in  phase-change  systems.  Consider  typical  solidification  systems.  Unless  the 
solidification  process  takes  place  under  microgravity  conditions,  say,  on  orbiting  space 
platforms,  gravity  is  always  present,  and  thermal  inhomogeneities,  which  also  are  always 
present,  give  rise  to  density  differences.  Density  gradients  also  can  be  caused  by 
compositional  inhomogeneities.  Depending  on  the  mutual  orientation  of  the  gravity  vector 
and  the  density  gradient,  buoyancy  driven  convection  can  be  set  up  which  can  significantly 
alter  the  temperature  and  composition  distribution  within  the  system  and  also  affect  the 
solid-liquid  interface  shape  and  dynamics.  If  there  is  a  free  surface  separating  the  liquid 
phase  and  a  surrounding  ambient  gas,  a  surface  tension  gradient  usually  will  exist,  caused 
by  thermal  and  composition  inhomogeneities,  which  will  exert  a  shear  stress  on  the  fluid 
phases,  leading  to  Marangoni  convection.  The  different  modes  of  convection  can  interact  to 
give  very  complex  flow  patterns  (Shyy  and  Chen  1993). 

Convection  also  can  be  caused  by  density  changes  during  phase  change.  Usually  the 
different  phases  have  quite  different  densities.  This  can  be  a  source  of  vigorous  convection 
as  the  fluid  moves  to  conserve  mass  (Shyy  1994).  Many  other  modes  of  convection  can  exist 
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and  interact,  such  as  forced  convection  (due  to  crystal  rotation,  etc.),  convection  due  to 
time-varying  electromagnetic  fields,  etc.  The  examples  shown  later  in  this  chapter  will  be 
confined  to  thermally  induced  buoyancy  driven  convection. 
4.5.4    Computational  Procedure 

The  solution  methodology  employs  a  finite  volume  formulation  in  the  context  of  the 
SIMPLE  algorithm  discussed  in  chapter  2.  Second-order  central  difference  schemes  are 
used  to  discretize  the  diffusion  and  the  convection  terms.  The  unsteady  terms  are  discretized 
using  the  implicit  first  order  Euler  differencing.  The  phase-change  source  terms  are 
discretized  in  a  manner  consistent  with  the  rest  of  the  terms  in  energy  equation.  Specifically, 

the  V  •  Vf  term  is  discretized  consistently  with  the  V  •  VT  term.  Further  details  of  the 

general  solution  procedure  can  be  found  in  Shyy  (1994). 

As  an  illustration  of  the  procedure,  the  melting  of  gallium  in  a  square  cavity  can  be 
considered.  For  this  case,  Gau  and  Viskanta  (1986)  have  reported  experimental  information 
of  interface  shapes  and  position  at  specified  time  instants.  The  shape  of  the  phase  change 
interface  is  determined  by  the  combined  effects  of  conduction  and  convection  driven  by 
buoyancy  forces.  It  will  be  shown  that  convective  effects  can  have  a  substantial  impact  on 
the  interface  shape. 

In  order  to  make  a  comparison  of  the  predicted  interface  shape,  at  various  time  instants, 
with  the  experimental  results  of  Gau  and  Viskanta  (1986),  the  calculations  were  carried  out 
with  the  following  parameters:  Ra  =  2.2  x  105  ,  Pr  =  0.0208  ,  St  =  0.042  .  A  value  of  E  = 
0.001  was  used  in  conjunction  with  the  T-based  method  to  ensure  numerical  stability.  The 
time  stepping  scheme  chosen  was  the  1st  order  backward  Euler  and  the  convection  scheme 
chosen  was  the  1st  order  upwind  method.  All  other  terms  were  discretized  by  2nd  order 
central  differences.  The  interface  location  was  determined  by  interpolation  for  the/  =  0.5 
contour.  Solutions  were  computed  on  a  21  x  21  uniform  grid  and  a  41  x  41  uniform  grid 
and  a  time  step  of  0.01  (1.43  seconds).  Numerical  results  have  been  shown  for  t  =  2  mins, 
5  mins,  6.667  mins  and  8  mins  in  Fig.  4.8  alongside  the  experimental  determinations  of  Gau 
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Figure  4.8.  Qualitative  comparison  of  predicted  interface  shapes  obtained  on  a 
41  X  41  uniform  grid  (solid  line)  with  experimental  data  (+)  of  Gau  and  Viskanta 
(1986)  for  the  melting  of  gallium  in  a  square  cavity. 

and  Viskanta,  after  translating  the  predicted  interface  to  account  for  the  undercooling  present 
in  the  experiment  (Lacroix  1989).  Reasonable  agreement  has  been  obtained  in  terms  of  the 
interface  shapes  at  the  various  time  instants.  The  effect  of  convection  is  to  transport  the  hot 
fluid  to  the  top  of  the  cavity  and  the  interface  moves  faster  near  the  top  of  the  cavity.  The 
shape  of  the  interface  is  distorted  by  the  convection  patterns  established  as  the  calculation 
progresses  in  time.  Thus,  it  is  evident  that  natural  convection  can  substantially  affect  the 
speed  and  the  shape  of  the  interface  from  the  planar  one  yielded  by  the  pure  conduction 
calculations. 

4.6  Effect  of  Gravity  Jitter 
Many  devices  in  spacecraft  contain  enclosures  filled  with  fluids.  In  the  following,  a 
problem  is  formulated  which  is  motivated  by  the  need  of  understanding  materials 
solidification  in  a  reduced  but  fluctuating  gravitational  field,  the  so-called  g-jitter  condition. 
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The  configuration  schematically  illustrated  in  Fig.  4.9  was  subjected  to  a  sinusoidal  variation 
in  the  magnitude  of  the  gravity  vector: 

g(t)   =  [1  +sin(o)t)]  (4.18) 

The  initial  condition  for  this  calculation  was  the  conduction  driven  solution  with  a  planar 
interface  at  the  x  =  0.5  location.  Calculations  were  then  conducted  with  an  angular 
frequency,  to  =  1  and  Rayleigh  number,  Ra  =  105,  based  on  the  maximum  value  of  the  gravity 
vector,  gmax-  A  uniform  grid  of  41  x  41  nodes  was  used  and  the  time  step  used  was  5% 
of  the  time  period  of  the  forcing  function. 

The  following  quantities  are  used  to  describe  and  assess  the  solution,  qj, ;  the  total  heat 
flux  through  the  hot  wall  (x  =  0),  qc  ;  the  total  heat  flux  through  the  cold  wall  (x  =  1).  A 
measure  of  the  instantaneous  interface  location  and  shape  was  obtained  by  sampling  the/= 
0.5  contour  at  a  total  of  5  points,  at  every  discrete  time  step.  The  coordinates  of  the  sampled 
points  may  be  described  as:  (a(t),0),  (b(t),0.25),  (c(t),0.5),  (d(t),0.75),  and  (e(t),l)  where 
a,  b,  c,  d  and  e  are  found  by  interpolating  for  the /=  0.5  contour,  along  y  =  0,  y  =  0.25,  y  = 
0.5,  y  =  0.75  and  y  =  1  respectively,  at  every  time  step.  A  schematic  describing  this  process 
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i — »      '  1 

y  =  0 

Figure  4.9   Schematic  showing  the  sampled  locations  for  the 
interface  position  at  each  time  instant  for  the  g-jitter  problem. 
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is  shown  in  Fig.  4.9. 

Figure  4. 10(a)  shows  the  time  history  of  the  heat  fluxes  qh  and  qc  along  with  their  phase 
relationships  to  the  forcing  frequency.  The  heat  flux  through  the  hot  wall,  qh,  is  mostly 
determined  by  convective  effects,  whereas  the  heat  flux  through  the  cold  wall,  qc,  is  affected 
by  conduction  only.  The  quick  response  of  convection  to  time  variations  in  the  body  force 
ensures  that  (a)  the  flux  qh  reaches  a  periodic  state  quickly  (about  2  time  periods),  and  (b) 
the  phase  difference  between  qh  and  the  forcing  function  is  negligible  for  the  chosen  forcing 
frequency,  co  =  1 .  The  heat  flux  qc  is  determined  by  conduction  effects  and  its  time  variation 
is  determined  by  the  change  in  boundary  conditions  caused  by  the  migration  of  the  T  =  0.5 
isotherm  that  constitutes  the  interface.  Hence,  qc  lags  the  forcing  frequency  by  90°.  The 
power  spectrum  of  qh,  and  qc  (10(b))  shows  the  existence  of  superharmonics,  but  the  energy 
content  of  the  superharmonics  are  at  least  two  orders  of  magnitude  less  than  that  of  the 
fundamental.  Thus  the  fundamental  frequency  is  the  dominant  one.  Figure  4.11(a)  shows 
the  time  history  of  the  interface  motion  at  the  sampled  points.  It  is  observed  that  the  portion 
of  the  interface  at  the  top  of  the  cavity  (y  =  1)  moves  towards  the  cold  wall  (x  =  1)  as  g 
increases  in  magnitude  and  convection  strength  increases,  whereas,  the  portion  of  the 
interface  near  the  bottom  of  the  cavity  (y  =  0)  moves  towards  the  cold  wall  (x  =  0).  The 
spectral  plots  (11(b))  show  that  the  fundamental  frequency  is  the  dominant  one  although 
higher  harmonics  exist. 

4.7  Float  Zone  Growth  of  NiAl 
In  the  following,  the  enthalpy  method  is  applied  to  the  float  zone  growth  of  NiAl  single 
crystals.  P-NiAl  is  an  intermetallic  that  is  currently  being  investigated  as  a  promising  high 
temperature  structural  material  for  application  in  the  next  generation  of  aircraft  engines  and 
structural  components.  NiAl  is  especially  attractive  because  of  its  low  density,  high  thermal 
conductivity,  high  melting  temperature,  superior  isothermal  and  cyclic  oxidation  resistance 
(Sen  and  Stefanescu  1991).  However,  it  has  two  major  drawbacks;  low  toughness  at  room 
temperatures  and  low  strength  at  high  temperatures  (Darolia  1 99 1 ).  These  deficiencies  need 
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Figure  4. 10.  (a)  Time  history  of  the  integrated  heat  fluxes  for  the  g-jitter  problem. 
The  solid  line  indicates  the  heat  flux  through  the  heated  wall  and  the  dashed  line 
indicates  the  heat  flux  through  the  cold  wall.  The  dotted  line  indicates  the  forcing 
function,  (b)  Spectral  plot  of  the  heat  flux  through  the  walls  of  the  cavity.  The  solid 
line  indicates  the  heat  flux  through  the  heated  wall  and  the  dashed  line  indicates  the 
heat  flux  through  the  cold  wall.  The  dotted  line  indicates  the  forcing  function. 
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(i)     Angular  Frequency  (ii)    Angular  Frequency 


Figure  4.11.  (a)  Time  history  of  the  sampled  interface  motion  at  every  time  step. 
The  letters  a,b,c,d,e  indicate  the  sampling  locations  as  shown  in  9  .  (b)  Spectral 
plots  of  the  sampled  interface  motion,  (i)  The  solid  line  refers  to  location  a  and 
the  dashed  line  to  location  e.  (ii)  The  solid  line  refers  to  location  b  and  the  dashed 
line  to  location  d.  In  both  cases,  the  dotted  line  refers  to  the  forcing  frequency. 
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to  be  resolved  before  NiAl  can  be  used  in  structural  applications.  The  current  study  is  part 
of  an  ongoing  effort  to  grow  high  purity  single  crystals  of  NiAl  with  controlled  compositions 
and  under  controlled  conditions.  This  would  be  the  starting  point  in  fundamental  and 
detailed  studies  directed  towards  understanding  the  behavior  of  this  material.  Figure  4.12 
shows  the  schematic  of  a  practical  float  zone  growth  system  for  producing  single  crystals 
of  NiAl.  Here  the  enthalpy  formulation  is  applied  to  a  practical  growth  system  as  shown  in 
the  schematic  in  Fig.  4. 12.  In  this  study  the  scope  is  restricted  to  the  convective  heat  transfer 
within  the  melt  and  its  impact  on  the  solid/liquid  interface  shape.  Therefore,  the  meniscus 
shape  and  the  heat  input  will  be  considered  given  quantities  during  the  course  of  the 
computations. 

4.7.1    Calculation  Procedure 

Issues  regarding  the  formation  and  stability  of  the  meniscus  shapes  have  been  discussed 
in  chapter  2.  Equation  (2.15)  with  associated  boundary  conditions  is  used  to  generate  the 
meniscus  shape.  In  the  current  set  of  calculations,  only  the  heat  transfer  within  the  melt  and 
its  effect  on  the  melt/crystal  and  melt/feed  interfaces  will  be  of  interest.  In  order  to  efficiently 
model  the  float  zone  with  adequate  resolution,  it  is  necessary  to  consider  only  half  the 
domain  and  apply  symmetry  conditions  at  the  centerline.  The  governing  equations  in 
nondimensional  form,  are 


(i)  continuity: 


3(Ru)  d(Rv) 

-^-  +  ^-  =  0  (4.19a) 


(ii)  radial  momentum: 
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Figure  4.12.  Illustration  of  the  float  zone  process.  The  float  zone  is  heated  by 
the  eddy  currents  induced  by  the  induction  coil  and  is  held  up  by  the  surface 
tension  between  the  melt  and  the  encapsulant.  The  quality  of  the  crystal  is 
controlled  by  the  thermofluid  transport  and  the  shape/velocity  of  the 
crystal/melt  interface. 

(iii)  axial  momentum 
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(iv)  energy: 
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and  the  boundary  conditions  are: 


106 


(i)  Top  and  bottom  boundaries: 

d2B 


dY2 

(ii)  Symmetry: 


=  0,  u  =  v  =  0  at  Y  =  0,  1.4  (4.20a) 


Q=f|=|f  =  0atR  =  0  (4-20b) 

(iii)  Free  surface  (melt  zone):  at  R  =  F(Y)  ,  0.2  <  Y  <  1.2 

(a)  Heat  flux:  q  =  V0  •  n  (4.20c) 

(b)  Marangoni  effect:  V(V  •  t)  ■  n  =  Ma  V0  •  t  (4.20d) 
(shear  stress  equals  the  surface  tension  gradient) 

(c)  Kinematic  condition:  V  •  n  =  0  (4.20e) 
(fluid  flow  is  tangential  to  the  free  surface) 

(iv)  Trijunction  points:  at  R  =  1,  Y  =  0.2  and  Y  =  1.2 

©  =  eiiquidus  =  0.2  (4.2CM) 

and  at  R  =  1,  Y  =  0  and  Y  =  1.4  0  =  0  (4.20g) 
where  the  variables  have  been  nondimensionalized  using  the  same  reference  scales  as  for 
equations  (3.15a  -  d).  The  treatment  of  the  additional  source  terms  arising  due  to  phase 
change  has  been  described  in  section  4.5.  The  geometry  and  boundary  conditions  have  been 
illustrated  in  Fig.  4.13. 

The  governing  equations  have  been  solved  using  the  pressure-correction  formulation 
in  nonorthogonal  body-fitted  coordinates.  The  convection  and  the  diffusion  terms  in  both 
the  momentum  and  the  energy  equations  have  been  discretized  using  second  order  central 
differences.  The  grid  points  were  clustered  towards  the  free  surface  and  for  the  computations 
involving  Marangoni  convection,  the  grid  points  also  were  clustered  towards  the  lower 
boundary  at  Y  =  0.  Only  steady  state  calculations  have  been  conducted.  Calculations  have 
been  carried  out  on  two  different  grids  -  a  coarse  grid  of  81  X  51  points  and  a  fine  grid  of 
161  X  101  points  -  to  assess  grid  independence  and  solution  accuracy.  The  fine  grid  system 
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Figure  4. 13.    Idealized  geometry  and  boundary  conditions. 


is  shown  in  Fig.  4. 14(a).  The  grid  has  been  clustered  towards  the  free  surface  to  obtain  better 
resolution  of  the  flow  and  heat  transfer  characteristics.  The  results  shown  in  the  following 
are  based  on  the  fine  grid  calculations. 
4.7.2    Results  and  Discussion 

It  will  be  shown  later  in  this  section  that  for  the  given  non-dimensional  parameter 
ranges,  surface  tension  gradients  generated  at  the  free  surface  will  play  a  dominant  role  in 
melt  convection  and  heat  transfer.  Based  on  the  thermophysical  properties  of  NiAl  and  the 
estimated  environmental  parameters  in  the  experimental  growth  configuration  of  Fig.  4.12, 
we  can  arrive  at  the  following  non-dimensional  parameters, 
(i)  Grashof  number,  Gr  =  2000 
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(ii)  Marangoni  number,   100  <  Ma  <  200  (dominant  source  of  convection) 

(iii)  Stefan  number,  St  =  0.1,  and 

(iv)  Prandtl  number,  Pr  =  0.1 

The  definitions  of  the  parameters  listed  above  are  based  on  those  discussed  in  section  2.3. 
The  heat  flux  is  assumed  to  be  input  normal  to  the  free  surface  and  is  estimated  based  on 
methods  outlined  in  chapter  2. 

4.7.3  Heat  Conduction 

A  calculation  was  carried  out  with  the  heat  transfer  taking  place  through  conduction 
effects  only.  The  purpose  was  to  carry  out  an  a  posteriori  verification  of  the  temperature 
scales  used  to  define  the  Grashof  and  the  Marangoni  numbers  in  subsequent  calculations. 
It  is  noted  that  due  to  the  curved  geometry  of  the  free  surface,  the  location  of  the  maximum 
temperature  shifts  towards  the  convex  portion  of  the  free  surface,  region  (B),  as  is  evident 
from  Fig  4. 14(b).  The  nondimensional  value  of  0max  is  2.47.  The  solid/liquid  interfaces 
show  moderate  curvature  at  the  axis  of  symmetry  (regions  C  and  D)  as  is  evident  from  Fig 
4.14(b). 

4.7.4  Thermocapillary  Convection 

Calculations  were  conducted  for  Gr  =  2000  and  two  Marangoni  numbers,  100  and  200. 
Figures  4. 15(a)  and  (b)  show  the  streamfunction  and  isotherms  respectively  for  the  case  of 
Gr  =  2000  and  Ma  =  100.  For  these  values  of  the  parameters,  the  surface  tension  gradient, 
in  the  form  of  Marangoni  convection,  plays  a  more  dominant  role  than  buoyancy,  in  the  form 
of  natural  convection,  in  determining  the  heat  transfer  characteristics.  Figures  4. 16(a)  and 
(b)  show  the  corresponding  quantities  for  Gr  =  2000  and  Ma  =  200.  It  is  evident  that  as  the 
Marangoni  number  increases,  the  solid/liquid  interfaces  show  increased  convexity  towards 
the  melt,  especially  near  the  axis  of  symmetry,  i.e.,  in  regions  (C)  and  (D).  A  clockwise 
convection  cell  is  formed  in  region  (B),  adjacent  to  the  convex  side  of  the  free  surface  and 
an  anticlockwise  cell  in  region  (A)  adjacent  to  the  concave  side  of  the  free  surface.  The 
location  of  the  maximum  temperature  is  pulled  downwards  towards  region  (B)  by  the 
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(a)      Grid:  161  X  101  points 


(b)     Conduction  solution 
Min:  0.003  Max:  2.472  20  Contours  at  uniform  intervals 


Figure  4.14  (a)  Illustrates  the  grid  used.  (A)  denotes  the  region  adjacent  to  the 
concave  side  of  the  free  surface;  (B)  denotes  region  adjacent  to  the  convex  side 
of  the  free  surface;  (C)  and  (D)  denote  the  regions  close  to  the  intersection  of  the 
solid/liquid  interfaces  and  the  centerline  as  shown. 

(b)  shows  the  conduction  solution.  The  solid/liquid  interfaces  are  mildly  convex 
towards  the  melt  in  regions  (C)  and  (D)  at  the  centerline  and  concave  close  to  the 
free  surface  at  regions  (A)  and  (B). 
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(a)    Ma  =  100,  Streamfunction 
Min:  -0.851  (-)  Max:  0.559  (+)  20  Contours 
-  denotes  clockwise  convection  and  +  anticlockwise  convection  cells. 


(b)    Ma  =  100,  Isotherms 
Min:  0.004  Max:  2.065  20  Contours 


Figure  4. 15  Ma  =  100.  (a)  Streamfunction:  The  bold  lines  indicate  the  solid/liquid 
interfaces.  The  increase  in  convexity  of  the  interfaces  in  regions  (C)  and  (D)  is 
evident.  There  is  a  clockwise  convection  cell  in  region  (B)  adjacent  to  the  convex 
side  of  the  free  surface  and  a  counterclockwise  cell  adjacent  to  the  concave  side  of 
the  interface  in  region  (A).  This  is  because  the  surface  tension  decreases  with 
temperature.  From  (b)  the  maximum  temperature  occurs  close  to  the  point  of 
inflexion  of  the  free  surface.  The  convex  region  of  the  free  surface  leads  to  a 
stronger  convection  cell. 


Ill 


(a)  Ma  =  200,  Streamfunction 
Min:  -1.349  (-)  Max:  0.711  (+)  20  Contours 
-  denotes  clockwise  convection  and  +  anticlockwise  convection  cells. 


(b)    Ma  =  200,  Isotherms 
Min:  0.006  Max:  1.887  20  Contours 


Figure  4.16  Ma  =  200.  The  convex  region  of  the  free  surface  leads  to  a  stronger 
convection  cell  in  region  (B).  Therefore  the  clockwise  cell  is  much  stronger  than 
the  counterclockwise  cell  in  region  (A).  The  location  of  maximum  temperature  is 
pulled  downwards  towards  (B)  by  the  clockwise  convection  cell.  The  magnitude 
of  maximum  temperature  has  decreased  significantly  due  to  increased  heat  transport 
by  the  convection  cells.  Decreased  heat  transfer  to  the  center  caused  increased 
convexity  of  the  solid/liquid  interfaces  near  the  centerline  and  increased  distortion 
near  the  free  surface  due  to  the  tangential  heat  transport  at  the  free  surface. 
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stronger  recirculation  of  the  clockwise  convection  cell,  which  in  turn  is  caused  by  the 
stronger  surface  temperature  gradients  there.  As  the  convection  strength  increases  with  the 
Marangoni  number,  the  magnitude  of  the  maximum  temperature  decreases  due  to  the 
increased  heat  transfer  effectiveness.  Increasing  the  strength  of  Marangoni  convection 
causes  the  heat  transport  to  be  tangential  to  the  free  surface.  As  a  consequence  of  the 
decreased  heat  transfer  in  the  normal  direction,  the  centerline  experiences  a  lower 
temperature.  The  increasing  temperature  inhomogeneity  due  to  this  type  of  convection 
pattern  results  in  increased  curvature  of  the  solid/liquid  interfaces  and  decreasing  volume 
of  the  melt.  The  trend  is  clearly  demonstrated  by  the  centerline  temperatures  shown  in  Fig. 
4.17(a),  which  shows  that  the  centerline  temperatures  decrease  with  increasing  Marangoni 
number,  and  from  Fig.  4.17(b),  which  shows  increasing  curvature  of  the  solid/liquid 
interfaces  as  a  result.  The  overall  effect  of  Marangoni  convection  is  to  increase  the 
temperature  inhomogeneity  of  the  melt,  especially  in  the  radial  direction,  causing  the 
solid/liquid  interfaces  to  experience  increased  curvature  towards  the  melt. 

The  above  series  of  calculations  demonstrate  the  application  of  the  enthalpy  method  and 
its  effectiveness  in  treating  problems  with  complex  geometry  and  involving  the  interaction 
of  the  solidification  interfaces  with  a  deformable  free  surface.  It  is  shown  that  surface  tension 
gradients  generate  strong  convection  in  the  vicinity  of  the  free  surface  which  significantly 
affects  the  heat  transfer  characteristics  and  interact  with  the  phase  change  dynamics.  The 
resulting  increases  in  the  curvature  of  the  crystal/melt  interface  will  have  a  significant  impact 
on  crystal  quality. 

4.8  Concluding  Remarks 
In  this  chapter,  alternative  numerical  formulations  of  the  enthalpy  method  have  been 
developed  and  assessed  for  the  computation  of  phase  change  of  pure  materials.  It  is  shown 
that  enthalpy  formulation  and  its  implementation  in  terms  of  the  T-based  method  can 
efficiently  capture  the  melt/solid  interface  with  adequate  resolution.  Convective  heat 
transfer  in  the  liquid  phase  has  also  been  taken  into  account.  A  variety  of  problems  relating 
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Figure  4.17.  Centerline  temperatures  are  reduced  in  the  melt  as  the  Marangoni 
number  increases  as  is  evident  from  (a).  The  solid-liquid  interface  is  drawn  inward 
towards  the  melt  at  the  centerline  but  pushed  towards  the  solid  near  the  free  surface 
as  the  Marangoni  number  increases,  as  shown  in  (b).  Marangoni  convection  causes 
the  heat  flux  to  move  tangentially  to  the  interface  and  decreases  the  heat  transported 
to  the  interior.  This  causes  the  interface  to  distort  and  become  convex  towards  the 
melt,  as  the  Marangoni  effect  becomes  stronger. 
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to  melting  and  solidification  have  been  considered  -  melting  of  gallium  in  a  cavity  and  float 
zone  growth  of  NiAl  single  crystals  -  all  of  which  are  of  immense  practical  importance.  The 
calculations  in  this  chapter  also  form  the  framework  for  a  more  extensive  and  detailed 
consideration  of  the  float  zone  process  in  the  next  chapter. 


CHAPTER  5 

MOVING  BOUNDARY  COMPUTATION  OF  FREE  SURFACE  FLOW 
AND  SOLIDIFICATION  IN  FLOAT  ZONE  PROCESSING 


5.1  Introduction 

In  this  chapter  a  numerical  methodology  has  been  developed  to  treat  fluid  flow  and  heat 
transfer  in  the  float  zone  process  for  single  crystal  growth.  A  moving  boundary  approach  has 
been  developed  to  handle  the  free  surface  between  the  melt  and  ambient  gas.  The 
solidification  interfaces  between  the  growing  crystal  and  its  melt  have  been  treated  using  the 
enthalpy  method  discussed  in  Chapter  4  as  opposed  to  the  explicit  interface  tracking 
employed  by  studies  described  earlier  in  Chapter  1.  The  enthalpy  method  makes  it 
unnecessary  to  explicitly  track  the  solid-liquid  interfaces,  thus  the  grid  generation 
requirements  are  less  stringent.  The  free  surface  has  been  computed  by  solving  the  dynamic 
form  of  the  Young-Laplace  equation  coupled  with  a  mass  conservation  condition.  The 
SIMPLE  algorithm  in  body-fitted  coordinates  (Shyy  1994)  has  been  utilized  to  solve  the 
convective  heat  transport  problem. 

5.2  Governing  Equations  and  Boundary  Conditions 
5.2.1    Transport  Processes 

The  Navier-Stokes  equations  of  continuity,  momentum  and  energy  can  be  written  in  a 

cylindrical  coordinate  system  and  specialized  to  the  axisymmetric  case  as  follows: 

.    .  d(pr)     d(pru)  3(prv) 

continuity:  JjjJ.  +  _i*L2  +  =  o  (5.1a) 

r-momentum: 
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at 


dr 


3y 


ap 
r  —  + 

9r 


dr 


2(xu 


Ar(u  -  ucast) 


(5.1b) 


y-momentum 


a(Qrv)  ^  a(Qruv)  ^  a(grvv) 
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energy: 
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(5.1c) 


(5. Id) 


where  the  source  terms  in  the  energy  equation  arise  from  the  enthalpy  formulation  for  phase 
change  problems  (Shyy  1994).  The  effect  of  property  variation  in  the  domain  has  been  taken 
into  account  and  the  properties  in  the  mushy  zone  are  given  by, 


thermal  conductivity: 


density: 


Lmush 


=  f  k,  +  (1  -  f)  ks 


Qmush  =  f  Qi  +  (1  -  f)  Qs 


(5.2a) 


(5.2b) 


where  the  subscripts  '1'  and 's'  indicate  the  liquid  and  solid  phases  respectively  and  f  is  the 
liquid  fraction.  The  D'arcy  law  term  is  used  to  drive  the  velocities  in  the  solid  to  equal  the 
casting  or  pulling  rate  given  by  vcast.  The  D'arcy  coefficient  is  given  by, 


A  =  C 


1  -  f 2 
f3  +  a 


(5.2c) 


where  C  is  a  large  number  and  q  is  a  small  number.  For  a  pure  material  they  can  be  chosen 
quite  arbitrarily  so  as  to  drive  the  velocity  in  the  solid  region  to  the  casting  rate  while 
retaining  numerical  stability.  The  liquid  fraction  is  modeled  as  a  linear  function  of  the 
temperature  within  the  phase  change  interval  and  is  set  to  zero  in  the  solid 
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f  = 


0  T  <  Ts 

T  -  T, 


2  £ 

1  ,     T  >  T 


Ts  <  T  <  T,  (5.2d) 


and  unity  in  the  liquid.  For  a  pure  material,  the  phase  change  interval,  2e,  is  chosen  so  as  to 
spread  the  numerical  mushy  zone  over  one  or  two  computational  cells  in  order  to  achieve 
numerical  stability.  As  regards  the  numerical  implementation,  a  choice  of  either  the  T-based 
or  the  H-based  method  is  available.  All  the  calculations  in  this  work  employ  the  T-based 
method. 

5.2.2    Boundary  Conditions 

5.2.2.1  Thermal  Boundary  Conditions 

Heat  transfer  to  the  system  is  controlled  by  radiation  and  convective  transfer  from  the 
ambient,  i.e., 

-  kVT  •  nlr=R  =  h(T  -  Tamb)  +  eo(T4  -  T^)  (5.3a) 

where  r  =  R  designates  the  surface  of  the  melt/crystal/feed  system.  Tamb  is  the  imposed 
temperature  distribution  due  to  the  heater.  Following  Duranceau  and  Brown  (1986)  and  Lan 
(1994),  Tamb  is  modeled  as  a  gaussian  distribution  given  by 


Tamb  =  (TP  "  Too)  expj^  -  i^-r2)  j  +  Too  (5.3b) 

where  Tp  is  the  heater  temperature  located  at  position  yp,  T<x>  is  the  ambient  temperature  far 
away  from  the  heater  and  a  is  the  half  width  of  the  distribution. 
5.2.2.2  Free  Surface  Treatment 

For  liquids  bounded  by  a  free  or  deformable  boundary  where  a  temperature  gradient 
exists  along  the  interface,  a  surface  shear  stress  is  generated  by  the  surface  tension  gradient. 
Thus,  the  force  balance  can  be  written  as, 

[i  V  (V  •  T)  •  n  =  |£  VT  •  T  (5.4a) 
where  the  surface  tension  gradient  is  proportional  to  the  temperature  gradient  and  gives  rise 
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to  the  normal  derivative  of  the  tangential  velocity  at  the  free  surface.  V  is  the  velocity  vector, 

n  is  the  unit  normal  vector  to  the  free  surface,  and  t  is  the  unit  tangent  vector  to  the  free 
surface.  The  convection  strength  is  described,  in  nondimensional  terms,  by  the  Marangoni 
number, 

dy 


Ma  = 


ATrb 

 b-  (5.4b) 


dT 

which  is  obtained  by  nondimensionalizing  Eq.  (5.4a)  with  a  temperature  scale,  AT,  length 
scale,  rb,  and  a  velocity  scale  given  by  a/r^,  where  a  is  the  thermal  diffusivity  of  the  melt 
and  |x  is  its  viscosity.  This  dimensionless  parameter  represents  the  strength  of  the  velocity 
gradients  generated  by  the  shear  stress  due  to  the  surface  tension  gradient. 

The  velocity  boundary  conditions  for  the  momentum  equations  are  obtained  by  solving 
Eq.  (5.4a),  in  conjunction  with  the  fact  that  the  meniscus  is  a  material  surface,  i.e., 

|£  +  u^  +  v|£  =  0  (5.5a) 
dt        dr        dy  v  ' 

where  F(x,  y,  t)  =  0  describes  the  equation  of  the  meniscus.  This  can  be  written  as, 

V  •  n  =  0  (5.5b) 
in  a  moving  grid  system  generated  to  fit  a  moving  and  deforming  meniscus.  Equations  (5.4a) 
and  (5.5b)  are  solved  simultaneously  to  obtain  the  boundary  velocities  at  the  free  surface. 
5.2.2.3  Symmetry 

Since  axisymmetry  is  assumed,  the  centerline  is  taken  to  a  symmetry  boundary.  Thus, 
at  r  =  0, 

u  =  0 
^  =  0 

dr     U  (5.6) 

f  =  0 
dr 
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5.2.3    Governing  Equations  for  a  Deformable  Meniscus 

The  meniscus  shape  obeys  the  Young-Laplace  equation  obtained  by  balancing  the 
hydrostatic  pressure,  hydrodynamic  pressure  and  viscous  normal  stress  terms  against  the 
surface  tension  force,  i.e., 


AQ  g  y  -  PL  -  Pd  +  2H 


dVn 

3n 


Y(T) 


1 


(1  +  r'2)3/2     r  (1  +  r'2)1/2 


(5.7a) 


where  r  =  r(y)  describes  the  free  surface  shape.  Here,  Aq  is  the  density  difference  between 
the  melt  and  the  ambient  fluid,  g  is  the  acceleration  due  to  gravity,  y  is  the  temperature 
dependent  surface  tension  between  the  melt  and  the  ambient  fluid.  The  term  on  the  right  hand 
side  represents  the  product  of  the  temperature  dependent  surface  tension  and  the  free  surface 
curvature.  The  left  hand  side  represents  the  hydrostatic  pressure,  the  hydrodynamic  pressure, 
Pd,  and  the  viscous  normal  stress  generated  by  the  fluid  motion.  In  addition,  there  is  a 
pressurization  parameter,  Pl,  that  must  be  determined  as  part  of  the  solution  either  by 
satisfying  a  contact  angle  condition  at  the  solidifying  interface  or  by  a  mass  conservation 
constraint. 

The  dynamic  form  of  the  Young-Laplace  equation,  Eq.  (5.7a),  can  be 
nondimensionalized  as  follows: 


BoY  -  WePd  -  PL  +  Ca 


dn  \ 


dy 

at\ 

dT 

R 


(1  +  R'2)3/2  R 


I  1 

(1  +  R'2)!/2j 


(5.7b) 

where  the  Bond  number,  Bo,  is  defined  as 

Bo=(Aq  g  r2)/(Yo)  (5.8a) 
describes  the  force  balance  between  the  hydrostatic  pressure  and  the  surface  tension,  the 
Weber  number,  We, 


We  =  (q  a2)/(Yo  rb) 


(5.8b) 


describes  the  relative  magnitudes  of  the  inertia  and  surface  tension,  and  the  Capillary 
number,  Ca, 
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Ca  =  (2  n  a)/(Yo  rb)  (5.8c) 
describes  the  relative  magnitudes  of  the  viscous  normal  stress  and  the  surface  tension.  Here, 
a  is  the  thermal  diffusivity  of  the  melt,  u,  is  the  viscosity  of  the  melt  and  f|,  is  a  characteristic 
length  scale  typically  taken  to  be  radius  of  the  melting  feed. 

Equation  (5.7a)  may  be  solved  with  the  following  boundary  conditions: 

(i)  Fixed  radius  at  the  melting  feed. 

r(y  =  hc)  =  rb  (5.9) 

and  either 

(ii)  Fixed  radius  at  the  solidifying  crystal, 

r(y  =  0)  =  rc  (5.10a) 
or  fixed  contact  angle  at  the  solidifying  crystal 

r'(y  =  0)  =  tan(ji/2  -  (J>c)  (5.10b) 
assuming  downward  pulling  as  is  the  case  in  the  problem  considered  subsequently.  Equation 
(10a)  has  been  used  in  all  the  cases  considered  in  this  study  for  the  purpose  of  comparison 
with  the  solutions  available  in  literature.  Here,  y  =  0  describes  the  location  of  the  melting 
feed  and  y  =  he  is  the  height  of  the  meniscus,  (j)c  is  the  contact  angle  between  the  meniscus 
and  the  solidifying  crystal  and  rc  is  the  radius  of  the  growing  crystal. 

5.3  Moving  Grid  Transformation 
Equations  (5.1a-d)  are  transformed  into  a  moving  grid  body  fitted  coordinate  system 
as  follows: 

The  transformation  may  be  expressed  as 

£  =  S(x,y,t) 

r\  =  n(x,y,t)  p.naj 

The  continuity  equation  in  nonorthogonal  coordinates  then  becomes, 

a  /T   x     d(orU)     d(prV)  n 

l(Jre)  +  "ir  +  "V  =  0  (5Jlb) 

where  r  is  the  radial  coordinate  which  is  identical  to  x  in  the  following  treatment.  The 
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separation  between  r  and  x  is  maintained  for  the  ease  of  application  to  a  two  dimensional 
Cartesian  coordinate  system  as  opposed  to  a  cylindrical  coordinate  system  considered  here. 


U  =  (u  -  x )  yn  -  (v  -  y )  Xt, 
V  =  (v  -  y )  xE  -  (u  -  x )  yE 


(5. lid) 


are  the  contravariant  velocity  components  and  x  and  y  are  the  radial  and  vertical 
components  of  the  grid  velocity  vector  respectively.  The  determination  of  the  Jacobian,  J, 
will  be  explained  subsequently.  The  momentum  equations  then  become 
r-momentum: 


d(Joru)  +  d(orUu)  d(orVu) 


at 


dt] 


_d 

all 


_a 


yr(q2ui  -  q3un) 


dp  dp 

^r(qlU^  -  q2uT)) 
+  Su(|,tl)  •  J  •  r 


(5.11e) 


y-momentum: 


3(Jrov)     a(grUv)  +  a(grVv) 


at 


an 


a(rp)  _  a(rp) 

xs  an     Xti  el 


+ 


_a_ 


^(qiv^-qiv,)]  +^[fr(q2vi-q3vT1) 
(O  -  Qref)  g  *  J  •  r  +  SV(|,T|)  '  J  •  r 


(5.1  If) 


along  x  and  y  directions  respectively,  where  x^,  y^  and  so  on  are  the  metrics  of  the 
transformation  and  qi,  q2  and  q3  are  defined  as 


qi  =  A  +  yl 

q2  =  x^  +  y^r, 

q3  =  A  +  yl 


(5.11g) 


Following  Shyy  (1994),  the  metric  terms  x^,  y^  and  so  on  are  evaluated  as  the  projected 
lengths  of  the  control  volume  in  the  physical  space  in  order  to  ensure  conservative  treatment. 
5.3.1    Determination  of  the  Jacobian 

The  Jacobian  of  the  transformation  is  defined  as, 
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J  =  x^t,  -  (5.12) 
which  can  be  interpreted  as  the  volume  of  the  computational  cell  (Shyy  1994).  However,  in 
a  moving  grid  system,  the  Jacobian  has  to  satisfy  the  geometric  conservation  law  which  can 
be  derived  by  considering  the  transformed  continuity  equation.  Consider  an  uniform 
flowfield  with  Q  =  1,  u  =  1  and  v  =  1 .  Such  a  flowfield  identically  satisfies  the  continuity 
equation  ( 1  a)  and  is  required  to  satisfy  the  transformed  continuity  equation  (5.11b).  Equation 
(5.11b)  then  becomes 


a[r(yx„   -xy,  )]  ^  jr  (xy§  -yx^  ) 


dt  d%  dT] 


=  0  (5.13) 


The  Jacobian  is  thus  computed  by  a  backward  Euler  difference  in  time  based  on  the 
geometric  quantities  such  as  the  radial  coordinate  and  the  metric  as  shown  above.  The 
methodology  is  an  extension  of  the  procedure  described  in  Shyy  et.  al.  (1996)  to  the 
cylindrical  coordinate  system. 

5.4  Computational  Methodology 

5.4.1    Meshing  Procedure 

For  moving  boundary  problems  like  the  float-zone  simulation,  it  is  essential  to  maintain 
an  adequate  resolution  in  the  liquid  zone,  even  as  the  grid  deforms  to  track  the  changing 
geometry.  Adequate  resolution  in  the  liquid  zone  is  necessary  in  order  to  (a)  stabilize  the 
calculation,  (b)  maximize  the  accuracy  of  the  convective  heat  transfer  and  phase  change 
computations,  and  (c)  accurately  capture  the  meniscus  deformation  and  faithfully  represent 
the  curved  geometry  of  the  domain.  At  the  same  time,  it  is  desirable  that  the  grid  generation 
process  be  fast  and  robust,  i.e.,  ensure  that  the  grid  distribution  is  smooth,  avoid  excessive 
skewness,  allow  good  control  of  the  grid  spacing  and  that  the  mapping  between  the  physical 
domain  and  the  computational  domain  is  invertible. 

In  this  work,  a  block  structured  approach  with  an  algebraic  grid  generation  procedure 
has  been  taken  to  enable  fast  grid  generation.  A  typical  float-zone  domain  can  be  divided 
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into  three  logical  subdomains:  (a)  the  crystal,  (b)  the  melt  and  (c)  the  feed  rod,  as  illustrated 
in  Fig.  5. 1 .  Subdomains  (a)  and  (c)  are  highly  stretched  to  fit  the  long  crystal  and  the  feed 
rod.  Thus,  the  material  in  subdomains  (a)  and  (b)  is  entirely  solid  during  the  solution  process. 
Subdomain  (b),  however,  encloses  the  melt  zone.  Thus,  phase  change  and  convective  heat 
transfer  occur  entirely  in  subdomain  (b).  It  is  also  to  be  noted  that  subdomain  (b)  has  an  0(  1) 
aspect  ratio  which  dictates  the  choice  of  the  interpolation  and  distribution  functions  chosen 
for  the  algebraic  grid  generation  procedure. 

Due  to  the  varying  and  sometimes  conflicting  requirements  listed  above,  the  grid  in  each 
block  or  subdomain  is  generated  separately  and  subsequently  assembled  together  into  a 
single  grid  system  to  continue  the  solution  procedure.  In  each  subdomain,  the  grid  is 
generated  algebraically  by  using  appropriate  interpolation  functions  and  distribution 
parameters,  checked  for  integrity  and  only  then  is  the  assembly  procedure  invoked  to  create 
a  single  grid  system.  The  algebraic  generation  procedure  is  relatively  inexpensive;  as 
opposed  to  the  solution  of  quasi-linear  elliptic  equations  and  has  been  found  to  be  robust  as 


Subdomain  (c):  the 

feed. 

Subdomain  (b):  the 

j 

melt. 

Subdomain  (a):  the 

growing  crystal 

Figure  5.1  Illustration  of  the  domain  and  its  partition  into  three  blocks 
for  the  purpose  of  grid  generation.  Block  2  enclosing  the  phase  change 
region  is  boxed  for  illustration.  The  extent  of  block  2  into  the  solid 
regions  can  be  adjusted  to  track  the  height  variation  of  the  zone. 
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well  in  the  present  problem.  However,  careful  choice  of  the  distribution  parameters  and 
interpolation  functions  is  necessary  to  ensure  smoothness  and  concentration  of  grid  lines 
close  to  the  interfaces. 

The  calculation  domain  is  divided  into  contiguous,  non-overlapping  blocks  and  the 
meshing  procedure  is  carried  out  separately  in  each  block.  For  convenience,  it  is  required 
that  each  block  consist  only  of  four  sides.  The  sides  may  be  specified  either  as  straight  lines 
or  by  sequence  of  (x,y)  pairs  describing  a  curve.  The  curve  may  be  specified  or  computed; 
the  only  requirement  is  that  it  does  not  intersect  itself.  For  instance,  in  the  melt  zone,  the 
meniscus  shape  is  computed  by  numerically  integrating  the  dynamic  form  of  the 
Young-Laplace  equation.  The  detailed  procedure  is  as  follows: 

(a)  Divide  the  domain  into  contiguous  blocks  of  four  sides  each.  Each  block  is  described  by 
a  non-intersecting  curve.  For  the  float-zone  calculations,  it  is  convenient  to  divide  the 
domain  into  three  blocks;  the  first  block  consisting  of  the  crystal,  the  second  block  consisting 
of  the  feed  rod  and  the  third  block  enclosing  the  melt  and  the  phase  change  region. 

(b)  For  each  side,  compute  the  arclength  from  the  sequence  of  (x,y)  pairs  describing  the 
curve.  The  required  number  of  grid  points  are  then  distributed  along  the  arclength  of  each 
side  as  illustrated  in  Fig.  5.2.  A  variety  of  choices  exist  for  the  distributing  function.  In  this 
study,  it  was  found  convenient  to  use  cubic  and  exponential  clustering  functions  depending 
on  the  geometry  of  the  individual  block.  The  cubic  polynomial  allows  a  fixed  number  of  grid 
points  to  be  distributed  with  specified  end  points  as  well  as  specified  end  spacings.  The 
exponential  clustering  is  suitable  for  highly  stretched  subdomains,  allowing  the  end  points 
to  be  matched  alongwith  the  spacing  at  one  end.  The  spacing  at  the  other  end  is  determined 
by  the  function.  The  exponential  clustering  is  described  by 

s  -  Sj      exp(a  ^)  -  1 


sn     si       exp(a)  -  1 


(5.14a) 


where  £  =  ^  _  *  ,  N  being  the  number  of  grid  points  that  need  to  be  accommodated  within 
Sj  <  s  <  sN.  The  coefficient,  a,  is  determined  by  specifying  the  spacing  at  one  end,  i.e.  as 
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Figure  5.2  Illustration  of  the  transformation  procedure  to  achieve  the 
desired  distribution  of  grid  points  along  the  arclength  of  the  boundary. 


the  root  of 


Asj        exp(aA^)  -  1 


(5.14b) 


sn  -  si       exp(a)  -  1 
where  the  root  is  determined  by  an  iterative  procedure.  The  quantities  |  and  T)  are  the 
transformed  and  normalized  coordinates  in  the  computational  plane  and  are  defined  as, 

i  -  1 


and 


i  = 


r\  = 


N;  -  1 

j  ~  1 
Nj-  1 


(514c) 
(14d) 


where  Nj  is  the  number  of  grid  points  along  the  £  direction  and  Nj  the  number  of  grid  points 

along  the  Tj  direction.  The  ability  to  match  the  spacing  at  one  end  allows  a  smooth  variation 
of  the  grid  spacing  between  the  different  subdomains. 

(c)  The  resulting  boundary  distribution  is  then  verified  by  ensuring  that  all  the  newly  formed 
arclength  elements  have  a  positive  length.  This  step  is  necessary  when  using  high  order 
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polynomials  such  as  the  cubic  clustering  because  of  potentially  oscillatory  behavior  between 
the  colocation  points. 

(d)  The  newly  formed  arclength  distribution  is  mapped  back  to  the  x-y  plane  via  cubic  spline 
interpolation.  This  completes  the  generation  of  the  boundary  distribution  for  each 
subdomain. 

(e)  The  internal  grid  distribution  is  generated  for  each  block  by  transfinite  interpolation.  For 

transfinite  interpolation,  use  of  linear  basis  functions  has  been  found  to  be  robust  and 

efficient.  The  quality  of  the  grid  is  determined  mostly  by  the  boundary  distribution  since 

using  low  order  interpolation  polynomials  does  not  give  rise  to  oscillatory  tendencies.  The 

disadvantage  is  that  it  is  does  not  provide  additional  control  of  the  interior  grid  distribution 

which  is  necessary  to  ensure  grid  orthogonality. 

7(^,  n)  =  ^  7(0,  n)  + (1-^)7(1,  n)  + 

T]  r(|,0)  +  (1  -  t|)  r(|,  1)  -  |  I)  7(0,0)  - 

|(1  -ti)?(0,1)-(1  -|)T|rU,0)-  (  } 

(1  -  I)  (1  -T|)  7(1,1) 

A  A 

where  7  =  x  i  +  y  j 

(f)  The  grid  is  checked  for  zero  Jacobians  and  convexity  of  the  control  volumes.  Loss  of 
convexity  implies  excessive  distortion  of  the  quadrilateral  that  constitutes  the  control 
volume.  Loss  of  convexity  is  tested  by  taking  the  cross  product  of  pairs  of  adjacent  sides  of 
the  control  volume  and  ensuring  that  the  angle  between  them  is  not  greater  than  n  radians. 
Appropriate  error  recovery  measures,  described  later,  are  taken  if  such  situations  are 
encountered. 

(g)  The  various  subdomains  are  combined  together  to  form  a  single  grid. 
5.4.2    Deformable  Meniscus 

The  meniscus  shape  is  computed  by  integrating  the  dynamic  Young-Laplace  equation 
described  earlier.  As  shown  in  Fig.  5.3,  the  meniscus  is  discretized  into  finite  domains  on 
the  underlying  grid  points.  The  dynamic  Young-Laplace  equation  is  then  integrated  over  the 
finite  domains  in  a  manner  analogous  to  the  discretization  of  the  transport  equations  over  the 
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Figure  5.3  Illustration  of  the  differencing  procedure  for  the  Young-Laplace 
equation. 

control  volume.  The  derivatives  in  the  equation  are  approximated  by  piecewise  linear 
interpolation  functions  effectively  giving  rise  to  a  second  order  central  difference  scheme 
for  all  the  terms  in  the  dynamic  Young-Laplace  equation. 

The  discretized  dynamic  form  of  equation  (5.7a)  can  now  be  written  in  the  form 


ap  Tj  =  as  rj_j  +  an  ri  +  1  +  b 
in  the  notation  popularized  by  Patankar  (1980),  where 

Y 


as  = 


Ayi-i/2(i  +  r'2) 


y  =  y 


j-1/2 


an 


Ayi  +  1/2(l+r'f/2 
ap  =  as  +  an 


y  =  yj+i/2 


and  the  source  term  is  given  by 


b  =  Qi  g 


yi±i  +  yL±  +  yi_y  f 

A  A  0  Jref 


-Pdi-PL  + 


Y 


rj  (1  +  r'2)1/2 

y=y: 


(5.16a) 


(5.16b) 


(5.16c) 


(5.16d) 


(5.16e) 


where  all  the  first  derivative  terms,  r',  are  evaluated  using  central  differences  over  the  nearest 
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half  grid  points  and  the  term  Pd .  includes  the  contribution  from  the  inertia  and  viscous  stress 

effects  due  to  melt  convection. 

The  equation  system  (5.16)  is  conveniently  solved  using  a  tridiagonal  solver.  However, 
since  the  coefficients  are  non-linear,  successive  substitution  is  employed  as  the  outer 
iteration.  During  the  course  of  the  iterations,  underrelaxation  is  employed  to  ensure 
convergence.  Underrelaxation  is  first  applied  to  the  non-linear  first  derivative  terms, 

r'k  =  (0r  5rk  +  (1  -  o)r)  r'k_1  (5.17a) 
where  k  denotes  the  iteration  number.  Underrelaxation  is  also  applied  to  the  shape  update 
as  follows: 

ap  =  ap/(om  (5.17b) 
b  =  b  +  (1  -  a)m)  ap  q  (5.17c) 
to  ensure  diagonal  dominance  before  solving  the  tridiagonal  system. 
5.4.2.1  Convergence  Criteria 

Convergence  criteria  for  terminating  the  non-linear  iterations  consists  of  examining  the 
residuals  summed  over  the  entire  field  as  well  as  the  1-norm  of  the  solution  variation  with 
iterations.  The  residual  is  always  scaled/normalized  by  the  solution  and  is  defined  as, 


Res  =  ^     an  rm  +  as  r^,  -  ap  r;  +  b  ^ 


(5.18) 


The  1-norm  is  defined  as, 


k-1 


£  |  rf  -  r| 

l|r|l="TT^T"  (5-19) 

i 

and  is  an  auxiliary  indicator  of  convergence.  Both,  the  residual  averaged  over  the  total 
number  of  control  volumes,  as  well  as  the  Lj-norm  are  required  to  be  less  than  10  ~4. 
5.4.2.2  Global  Mass  Conservation 

Meniscus  shapes  are  determined  by  adjusting  the  value  of  PL  so  as  to  satisfy  the 
constraint  of  global  mass  conservation.  In  this  study,  a  target  domain  mass  is  determined  by 
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considering  the  domain  to  be  initially  a  solid  cylinder.  As  melting  progresses,  a  mass  deficit 
builds  up  due  to  density  differences  between  the  melt  and  the  solid.  The  meniscus  shape  is 
adjusted  to  reduce  the  mass  deficit  to  zero  via  manipulation  of  PL. 
First,  the  mass  deficit  is  determined, 


AM  =  Mdomain  -  [    q  dV  (5.20) 

J  Vol 


The  meniscus  is  then  required  to  satisfy  a  volume  constraint, 


•I 


dy  +  ^  (5.21) 


The  required  value  of  PL  is  determined  by  finding  the  root  of 

^(Pl)  =  Vm  -  jt  I  r  dy  =  0  (5.22) 


5.4.2.3  Root  Finding  Procedure 

(i)  The  root  is  bracketed  between  two  extreme  values  of  PL  by  a  linear  search  for  a  zero 
crossing  of  ^(Pl).  It  is  helpful  to  carry  out  exploratory  computations  to  arrive  at  reasonable 
guess  values  for  the  step  size  as  well  as  the  initial  guess  for  PL. 

(ii)  The  secant  method  is  used  to  locate  the  root.  If  a  paralogical  situation  is  encountered 
where  the  secant  method  does  not  converge,  the  procedure  automatically  switches  over  to 
the  bisection  method. 

5.4.2.4  Remeshing  Procedure 

Once  a  meniscus  profile  is  obtained  satisfying  the  dynamic  form  of  the  Young-Laplace 
equation  and  the  global  mass  conservation  constraint,  the  domain  has  to  be  remeshed  as  the 
grid  moves  to  fit  the  meniscus.  It  is  possible  to  underrelax  the  grid  movement  in  order  to 
stabilize  the  transport  calculations.  Again,  an  algebraic  procedure  is  adopted  whereby  the 
x-coordinates  of  the  control  volumes  are  determined  by  using  transfinite  interpolation, 
(i)  Carry  out  linear  transfinite  interpolation  to  determine  the  x-coordinates  of  the  grid 
points  using  the  x-component  of  Eq.  (15). 


130 


(ii)  Check  the  grid  for  zero  Jacobians  and  loss  of  convexity  of  the  control  volumes.  The  error 
recovery  procedure  is  described  next. 

5.4.3  Error  Recovery 

It  is  possible  to  encounter  zero  Jacobians  during  the  iterative  process  if  the  height  of  the 
zone  decreases  so  that  the  grid  gets  excessively  squashed.  In  such  a  case,  the  algorithm  prints 
a  warning  message  and  continues  the  calculations  without  updating  the  grid.  To  prevent  this 
from  occurring  it  is  possible  to  specify  a  minimum  zone  height  below  which  the  vertical 
component  of  grid  motion  is  suppressed.  Loss  of  convexity  usually  occurs  if  the  meniscus 
computation  diverges  due  lack  of  existence  of  solutions  in  the  required  parameter  range  - 
especially  if  the  zone  height  becomes  excessively  large.  An  appropriate  parameter  range  can 
be  approximately  determined  by  solving  the  static  form  of  the  Young-Laplace  equation  to 
generate  an  existence  map  over  the  desired  parameter  range  and  utilizing  such  a  map  to 
generate  appropriate  initial  guesses.  Such  a  procedure  will  be  described  in  detail  in  the  results 
section. 

5.4.4  Grid  Underrelaxation 

It  will  be  demostrated  later  that  grid  underrelaxation  is  necessary  in  order  to  ensure  a 
smooth  and  oscillation  free  convergence  history.  The  vertical  and  radial  components  of  the 
grid  deformation  are  underrelaxed  separately  and  the  corresponding  underrelaxation 
parameters  are  denoted  by  cor  and  0)y  respectively.  Typical  values  of  cor  range  from  0.4  to  0.7 
depending  on  the  estimated  height  of  the  melt  zone.  The  vertical  component,  however,  needs 
to  be  underrelaxed  substantially.  The  value  of  ooy  has  been  varied  from  0  (no  deformation 
in  the  y-direction)  to  1  (no  underrelaxation)  in  an  attempt  to  quantify  its  effect  on  the  solution 
history.  Typical  values  used  are  in  the  range  of  0.05  to  0. 1  to  ensure  convergence  of  the 
solution.  Care  is  taken  to  ensure  that  the  melt  zone  is  completely  enclosed  by  the  refined  grid 
region  comprising  subdomain  (b)  as  illustrated  in  Fig.  1  by  tracking  the  location  of  the 
subdomain  boundaries  with  respect  to  the  phase  boundaries  during  the  iterative  process  and 
also  by  an  a  posteriori  estimation  after  the  solution  has  been  obtained. 
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5.5  Results  and  Discussion 

5.5.1  Test  Problem 

The  test  problem  used  to  develop  and  test  the  algorithm  is  the  float  zone  growth  of 
Sodium  Nitrate  crystal  (NaNC^)  crystals  heated  using  a  ring  heater.  The  steady  state 
solutions  for  pure  heat  conduction,  meniscus  shapes,  natural  and  Marangoni  convection  are 
available  in  literature  (Lan  and  Kou  1990a,  1991a)  and  Lan  (1994).  These  solutions  provide 
a  reference  point  for  comparison  and  assessment. 

An  initially  cylindrical  rod  of  NaNC>3  is  heated  at  the  center  by  radiation  from  a  ring 
heater  which  causes  it  to  start  melting  and  form  a  float  zone.  The  meniscus  deformation 
depends  on  the  length  of  the  molten  zone  and  also  on  the  density  change  during  melting  as 
the  meniscus  deforms  to  conserve  mass.  The  rod  is  then  pulled  downwards  to  simulate  the 
crystal  growth  processes.  Since  the  length  of  the  rod  is  set  to  40  times  its  radius,  the  melt  zone 
occupies  only  a  small  fraction  of  the  domain.  The  physical  properties  of  NaNC>3  are  shown 
in  table  5.1.  Other  relevant  geometrical  and  physical  parameters  are  presented  in  table  5.2. 
In  order  to  accommodate  the  enthalpy  method,  the  melting  temperature  is  distributed  over 
a  finite  phase  change  interval  given  by 

2e  =  0.2  K  (5.23) 

5.5.2  Nondimensional  Parameters 

The  material  properties  and  the  geometrical  parameters  presented  in  tables  5. 1  and  5.2 
above  enable  the  estimation  of  important  nondimensional  parameters  for  the  system.  Based 
on  an  estimated  temperature  scale,  AT  =  Tmax  -  Tm  =  20K,  length  scale,  rb  =  0.2cm, 
and  the  material  properties  given  in  table  5. 1 ,  the  following  nondimensional  parameters  can 
be  obtained: 

Prandtl  number  Pr  =  a   =  9.2 

g(3ATr3 

Grashof  number   — -  =  411 

v2 

Rayleigh  number  Ra  =  Gr  Pr  =  378 1 
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Table  5.1.  Physical  properties  of  NaNC>3 

Melting  temperature,  Tm 

579.8  K 

Latent  heat,  L 

182  J g-1 

Thermal  conductivity  (solid), 

5.65  X  10~3  +  33.5  (T  -  503)  X  10"7  W  cm"1  Kr1 

Thprmal  ronHnctivitv  Hinnid^ 

1  11^-1  11  iLll    \*/ UUUUV/ LI  V  1 L  V  IIIUUIUL 

5  65  X  1 0-3  +  33  5  <T  -  T     V  1  Cr1  W  nrr1  K~] 

Specific  heat,  Cps  =  Cpi 

1.255  +  2.18  (T  -373)  X  10~3  J  g"1  K"1 

Surfarp  tpnsion  praHipnt  r)v/r)T 

-  0  056  dvn  cm-1  K_1 

uy  ii  win  iv 

OllllClV^V^   l^ll.llv'll,  y 

1 19  96  -  0  056  f  T  -  T   )  dvn  rm-1 

uyiiaiiuc  VlaLUSliy,  |X 

u.ujuz  —  i.jjj  x  io    \i  —  i g  cm  s 

Thermal  expansion  coefficient, 
P 

6.6  X  10"4  K-1 

Density  (solid),  qs 

2.118  gem-1 

Density  (liquid),  Q\ 

1.904  gem-1 

Table  5.2.  Geometrical  and  physical  parameters 

rb  (radius) 

0.2  cm 

yp  (heater  location) 

4.0  cm 

a  (width  of  heater  profile) 

0.3  cm 

Toe 

298  K 

Tp  (heater  temperature) 

813  K,  913  K 

e  (emissivity) 

0.7 

h  (heat  transfer  coefficient) 

4.6  x  10^  W  cm-2  K-1 

Marangoni  number 
Stefan  number 


Bond  number 
Capillary  number 
Weber  number 


Ma  =  4752 
Cn  AT 

St  = 


AH 


0.3 


AQgr2, 
Bo  =  -rr-±   =  0.61 


Ca 


Yo 


3.8  x  10"6 


We  = 


qccz 


=   =  2.4  x  10"7 
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Biot  number  =  0.02 

k 

earhAT3  , 

Radiation  number   r         =  3.7  x  10  ~6 

k 

5.5.3    Steady-State  Calculations 

A  conduction-capillarity  solution  obtained  by  setting  the  heater  temperature  to  913  K 
has  been  chosen  to  present  some  of  the  features  of  the  algorithm  and  its  implementation,  and 
to  characterize  the  interaction  between  the  grid  movement,  meniscus  deformation  and  heat 
conduction  terms  in  the  discretized  energy  equation. 
5.5.3.1  Estimation  of  Meniscus  Shape  and  P]_ 

The  meniscus  calculations  need  a  starting  procedure  to  determine  the  initial  guesses  for 
the  melt  zone  height  and  the  corresponding  guess  value  of  PL.  The  procedure  is  as  follows: 

(i)  Determine  an  initial  estimate  of  the  melt  zone  height.  This  was  done  by  computing  a  pure 
heat  conduction  solution  for  a  cylindrical  geometry,  i.e.,  with  the  meniscus  shape  held  fixed. 

(ii)  Using  the  guessed  height,  the  static  Young-Laplace  equation  is  solved  using  an  adaptive 
fifth  order  Runge-Kutta  method  and  the  shooting  technique  to  generate  all  possible  profiles 
over  the  entire  range  of  PL  for  which  solutions  exist.  The  relative  local  error  (Shyy  1994) 
is  maintained  below  10-6  by  comparing  with  the  corresponding  fourth  order  Runge-Kutta 
technique  and  varying  the  step  size  accordingly. 

(iii)  Figure  5.4(a)  shows  the  map  of  melt  volume,  Vm,  versus  PL  over  the  entire  range  of 
available  solutions.  It  is  evident  that  over  the  range  of  interest,  two  solutions  exist  even  for 
the  static  Young-Laplace  equation.  Figure  5.4(b)  shows  the  meniscus  profiles  satisfying  the 
volume  constraint  with  a  relative  tolerance  of  10~2. 

(iv)  Figures  5.4(a)  and  (b)  are  used  to  guide  the  initial  guesses  for  Pl  and  the  limiting  values 
for  the  meniscus  volume  for  the  float  zone  computations.  Such  a  procedure  ensures  operation 
in  the  range  where  the  solutions  to  the  Young-Laplace  equation  exist. 
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Figure  5.4.  Determination  of  initial  guesses  for  melt  zone  height  and 
Pl  using  the  shooting  method  with  Runge-Kutta  fifth  order  integration. 

(a)  Profiles  satisfying  the  static  Young-Laplace  equation  and  the  given 
melt  volume  with  a  relative  tolerance  of  10-2  and  local  error  less  than 

10-6. 

(b)  A  map  of  melt  volume  versus  Pl  over  the  entire  range  of  Pl  for  which 
the  static  Young-Laplace  equation  can  be  integrated.  It  is  seen  that  over 
a  significant  range  of  Pl,  two  solutions  exist  satisfying  a  given  melt 
volume  constraint. 
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5.5.3.2  Conduction-Capillarity  Solution 

Two  sets  of  calculations  with  different  initial  guesses  and  grid  underrelaxation 
parameters  were  carried  out  to  investigate  the  coupling  between  the  grid  movement, 
boundary  deformation  and  heat  conduction.  First,  a  conduction  solution  was  obtained  with 
the  heater  temperature  set  to  913  K  and  by  imposing  a  cylindrical  shape  on  the  meniscus. 
This  was  then  used  as  an  initial  guess  for  the  conduction-capillarity  problem  which  allowed 
the  meniscus  shape  to  deform  in  accordance  with  the  Young-Laplace  equation.  No 
underrelaxation  was  used  when  updating  the  meniscus  shape  or  the  grid.  Figure  NO  TAG 
shows  the  convergence  characteristics  of  this  procedure.  Figure  5.5(a)  shows  the  history  of 
the  monitored  quantities,  Hm  -  the  zone  height  at  the  free  surface,  and  Hc  -  the  zone  height 
at  the  centerline.  Figures  5.5(b)  and  (c)  show  the  history  of  Pl  and  the  melt  volume,  Vm, 
respectively.  The  residual  of  the  energy  equation  and  the  1 -norms  of  the  temperature,  x,  and 
y  components  of  the  grid  deformation  are  shown  in  Fig.  5.5(d-f).  The  highly  oscillatory 
trend  in  the  solution  characteristics  is  evident  and  this  procedure  cannot  yield  a  solution. 

The  second  set  of  calculations  used  an  initially  solid  rod  at  room  temperature,  300  K, 
as  the  initial  guess.  Figure  5.6  shows  the  effect  of  various  grid  underrelaxation  factors  on  the 
convergence  characteristics.  Figures  5.6(a-c)  shows  the  history  of  Hm  for  grid 
underrelaxation  factors  over  a  range,  0  <  coy  <  1.  0)y  =  1  indicates  that  the  vertical 
component  of  grid  deformation  was  not  underrelaxed  and  toy  =  0  indicates  that  the  vertical 
component  of  grid  deformation  was  completely  suppressed.  Figures  5.6(d-f)  show  the 
corresponding  trends  for  the  quantity  Hc.  It  is  evident  that  grid  underrelaxation  in  the  range, 
0. 1  <  (Dy  <  0.05  is  sufficient  to  suppress  the  oscillatory  behavior.  Complete  suppression 
of  the  vertical  component  of  grid  deformation  yields  the  smoothest  but  the  slowest 
convergence  characteristics. 

Figure  5.7  shows  the  solution  features.  It  is  observed  that  the  meniscus  profile  and  the 
zone  height  are  very  close  to  the  profile  shown  in  Fig.  5.4a.  The  hydrostatic  pressure  causes 
the  melt  to  bulge  outwards  near  its  base  and  to  'neck'  inwards  near  the  top.  Figure  5.7(a) 
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Convergence  history  for  conduction  case  with  Tp  =  640°C 
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Figure  5.5.  Convergence  history  for  the  conduction  case  with  Tp 
=  9 13  K.  Effect  of  initial  guess  and  underrelaxation  parameters  on 
the  convergence  history.  A  conduction  solution  was  obtained  by 
imposing  a  cylindrical  meniscus  shape.  The  converged  solution 
was  then  used  as  the  initial  guess  for  the  conduction  -  capillary 
solution.  The  grid  motion  and  boundary  deformation  were  not 
underrelaxed,  i.e.,  toy  =  wr  =  1. 

shows  the  computed  grid  and  the  free  surface  shape.  Figure  5.7(b)  shows  the  corresponding 
isotherms.  The  phase  interfaces  are  illustrated  by  the  thick  lines  in  Fig.  5.7(b).  The 
melt/crystal  interface  at  the  base  of  the  zone  is  quite  flat  whereas  the  melt/feed  interface  at 
the  top  is  slightly  concave  towards  the  melt.  Figure  5.8  shows  the  temperature  profiles  along 
the  entire  height  of  the  domain.  The  temperature  profile  at  the  surface  is  shown  by  the  solid 
line  and  the  profile  at  the  centerline  is  shown  by  the  dashed  line.  It  is  observed  that  the 
temperature  distribution  in  the  solid  regions  becomes  essentially  one  dimensional  away  from 
the  melt  zone.  This  confirms  the  assumption  made  by  Duranceau  and  Brown  (1986)  in  their 
conduction-capillarity  model  which  imposed  a  one  dimensional  temperature  distribution  in 
regions  far  from  the  melt  zone  in  order  to  reduce  the  computational  load.  The  present 
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Figure  5.6.  Convergence  history  for  the  conduction  case  with  Tp  =  913  K. 
Effect  of  initial  guess  and  underrelaxation  parameters  on  the  convergence 
history. 
cor  =  0.7. 

(a)  -  (c)  History  for  Hm  for  various  ooy. 

(d)  -  (e)  History  for  Hc  for  the  corresponding  toy. 

In  (d)  and  (e),  u)y  =  0  indicates  that  only  the  x  component  of  grid  deformation 
is  permitted.  The  fixed  grid  solution  is  slower  to  converge,  but  the  solution 
history  is  smoother.  The  initial  guess  was  a  solid  cylindrical  rod  at  300  K. 
Calculations  on  grid  of  41  X  101  points. 


computations  help  to  justify  this  assumption  and  also  designate  the  regions  in  which  a  one 
dimensional  distribution  can  be  prescribed. 

5.6  Grid  Refinement  Study  and  Assessment 
A  series  of  calculations  have  been  conducted  with  parameters  chosen  so  as  to  facilitate 
comparison  with  the  solutions  of  Lan  (1994).  The  heater  temperature  was  set  to  813  K. 
Starting  with  the  conduction-capillarity  solution,  the  effects  of  buoyancy  induced 
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Figure  5.7.  Heat  conduction  solutions  with  Tp  =  913  K. 

(a)  Grid  and  boundary  shape. 

(b)  Isotherm  distribution  showing  the  isotherm  distribution  on  the  grid  of 
41  X  101  points.  The  interfaces  are  denoted  by  the  thick  lines. 


convection  and  thermocapillary  convection  have  been  successively  incorporated.  For  all 
computations,  the  convergence  history  of  important  quantities  such  as  the  height  of  the  zone 
at  the  free  surface,  Hm,  the  height  of  the  zone  at  the  centerline,  Hc,  and  the  heat  flux  balance 
has  been  monitored.  These  calculations  employ  the  second  order  accurate  central  difference 
scheme  for  the  convection  and  diffusion  terms  in  the  discretized  form  of  the  Navier-Stokes 
equations.  These  calculations  were  carried  out  to  assess  the  present  computations  with  those 
of  Lan  (1994)  in  terms  of  accuracy  and  numerical  diffusion.  The  study  by  Lan  (1994)  and 
the  other  works  cited,  have  used  a  first  order  upwind  scheme  to  discretize  the  convection 
terms. 
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Temperature  profiles  for  conduction  solution  with  Tp  =  913  K 
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Figure  5.8.  Temperature  profiles  along  the  height  of  the  domain  at  the  surface 
and  at  the  centerline.  The  profiles  show  that  a  one  dimensional  temperature 
distribution  is  established  in  the  feed  and  crystal  regions  away  from  the  melt 
region.  This  confirms  the  one  dimensional  approximation  used  by  Duranceau 
and  Brown  (1986)  in  the  regions  far  from  the  melt. 

5.6.1  Grid 

Two  grids  have  been  used  for  the  accuracy  assessment,  a  grid  of  41  X  101  points  and 
a  refined  grid  of  81  X  181  points. 

(T)41  x  101  grid:  This  grid  has  41  points  in  the  radial  direction  with  a  spacing  of  0.002  cm 
at  the  meniscus  and  0.008  cm  at  the  symmetry  boundary  distributed  by  fitting  a  cubic 
polynomial.  In  the  vertical  direction,  the  phase  change  zone  has  61  points  distributed  by  a 
cubic  polynomial  with  the  spacings  at  each  end  set  to  0.007  cm.  The  feed  and  crystal  zones 
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have  21  points  each,  distributed  according  to  the  exponential  clustering  by  matching  the 
spacings  to  the  phase  change  zone.  The  grid  in  the  phase  change  zone  extends  beyond  the 
location  of  the  trijunction  points  by  0.02  cm. 

(ii)  81  x  181  grid:  This  grid  has  81  points  in  the  radial  direction  with  a  spacing  of  0.001 
cm  at  the  meniscus  and  0.004  cm  at  the  symmetry  boundary  distributed  by  fitting  a  cubic 
polynomial.  In  the  vertical  direction,  the  phase  change  zone  has  121  points  distributed  by 
a  cubic  polynomial  with  the  spacings  at  each  end  set  to  0.003  cm.  The  feed  and  crystal  zones 
have  31  points  each,  distributed  according  to  the  exponential  clustering  by  matching  the 
spacings  to  the  phase  change  zone.  The  grid  in  the  phase  change  zone  extends  beyond  the 
location  of  the  trijunction  points  by  0.03  cm. 
5.6.1.1  Heat  Conduction  Calculations 

The  initial  guess  was  a  solid  cylindrical  rod  at  300  K.  Initially,  the  meniscus  deformation 
was  not  permitted  and  a  cylindrical  meniscus  shape  was  imposed.  A  converged  solution  was 
obtained  for  this  configuration  when  the  zone  height  stopped  varying  and  the  residual  for 
the  energy  equation,  averaged  over  the  control  volumes,  was  less  than  10-5.  This  was  used 
as  the  initial  guess  for  the  subsequent  set  of  calculations  with  the  meniscus  shape  computed 
by  integrating  the  Young-Laplace  equation.  The  zone  height  is  a  very  sensitive  indicator  of 
solution  convergence  since  it  occupies  a  small  fraction  of  the  overall  length  of  the  domain. 
Calculations  have  been  carried  out  with  the  meniscus  shapes  computed  under  both  normal 
and  under  microgravity  (Og)  conditions  to  assess  the  impact  of  meniscus  shape  on  the 
solution.  Figure  5.9  shows  the  convergence  history  of  the  solution  for  the  normal  gravity 
case.  Figure  5.9(a)  and  (b)  show  the  history  of  the  zone  height  -  Hm  (at  the  surface)  and  Hc 
(at  the  centerline)  for  the  grid  of  41  X  101  points  and  for  the  grid  of  81  X  181  points, 
respectively.  Figure  5.9(c)  and  (d)  show  the  history  of  the  heat  flux  balance  in  the  domain. 
These  plots  establish  the  convergence  of  the  solution. 

Figure  5.10  (i)  and  (ii)  show  the  computed  shape  and  converged  grid  distribution 
obtained  on  the  4 1  x  101  grid  and  the  8 1  X  181  grid  respectively.  The  zone  heights  and  the 
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Figure  5.9  Convergence  history  for  the  conduction  solution  at  lg.  (a)  and 
(b)  height  of  the  zone  at  the  surface  and  at  the  centerline  for  the  grid  of  41 
X  101  points  and  81  X  181  points  respectively,  (c)  and  (d)  balance  of  the 
incoming  and  outgoing  heat  fluxes  for  the  grid  of  41  X  101  points  and  81 
X  181  points  respectively. 


maximum  temperature,  Tmax,  agree  to  within  0.6%  and  0.2%  respectively,  establishing  the 
grid  independent  nature  of  the  converged  solutions.  Figure  5.11  shows  the  corresponding 
isotherms  and  the  solid  -  liquid  interfaces  are  shown  by  the  thick  lines.  The  menicus  has  a 
pronounced  bulge  at  the  base  of  the  zone  due  to  the  combined  action  of  gravity  and  the 
density  change  due  to  phase  change.  The  solid  liquid  interfaces  show  mild  convexity  towards 
the  melt.  The  solution  on  the  41  x  101  grid  agrees  with  the  solution  of  Lan  (1994),  the 
difference  in  Tmax  being  less  than  0.6%.  Figure  5.12  (a)  and  (b)  show  that  the  temperature 
profiles  obtained  on  the  41  x  101  grid  and  the  81  X  181  grid  are  identical.  Figure  5.12(c) 
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Hm  =  0.328  cm,                                       Hm  =  0.328  cm, 
(j)   H,  =  0.236  cm.  (")  Hr  =  0.234  cm.  

Figure  5.10  Computed  grid  and  boundary  shape  for  the  conduction  solution  at  lg.  (i) 
Grid  of  41  X  101  points,  (ii)  Refined  grid  of  81  X  181  points.  Only  the  phase  change 
region  and  its  immediate  vicinity  is  shown  since  the  actual  height  of  the  rod  is  40  times 
the  radius.  The  zone  heights  agree  to  within  0.2%  at  the  surface  and  0.6%  at  the 
centerline. 


.59  K, 


Tmin  =  298.00  K,  PL  =  994.79  Tmin  =  298.00  K,  PL  =  988.37 

dynes/cm2  dynes/cm2 

Figure  5.11  Isotherms  for  the  conduction  solution  at  lg.  (i)  Grid  of  41  X  101  points 
and  (ii)  Refined  grid  of  81  X  181  points.  The  Tmax  agree  to  within  0.02%.  The  41 
X  101  grid  solution  differs  with  Lan  (1994)  by  0.6%. 
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Temperature  profiles  for  conduction  solution 
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Figure  5.12  (a)  and  (b)  Temperature  profiles  showing  comparison 
between  coarse  grid  of  41  X  101  points  and  fine  grid  of  81  X  181  points 
along  the  entire  height  of  the  domain,  (a)  Temperature  profile  along  the 
surface,  (b)  Temperature  profile  along  the  centerline. 
(c)  and  (d)  Comparison  of  Temperatures  along  the  surface  and  the 
centerline.  (c)  41  X  101  grid,  (d)  81  X  181  grid. 


and  (d)  show  that  the  temperature  distribution  in  the  feed  and  crystal  regions  becomes 
one-dimensional  far  away  from  the  zone. 

Figures  5.13-5.16  correspond  to  the  microgravity  condition.  The  most  significant 
impact  of  microgravity  on  the  conduction  solution  is  on  the  meniscus  shape.  The  meniscus 
shape  now  resembles  a  circular  arc  due  to  the  absence  of  gravity  and  the  solution  displays 
symmetry  about  the  heater  location.  Figure  5.13  establishes  the  convergence  of  the  solutions 
on  the  41  x  101  grid  and  the  81  x  181  grid.  Figure  5.14  and  5.15  show  the  computed 
boundary  shape,  grid  distribution  and  the  isotherms.  It  is  worth  noting  that  the  value  of  Pl 
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Solution  history  for  conduction  solution  with  deforming  meniscus  at  Og 
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Figure  5.13  Convergence  history  for  the  conduction  solution  at  Og.  (a) 
and  (b)  height  of  the  zone  at  the  surface  and  at  the  centerline  for  the  grid 
of  41  X  101  points  and  81  X  181  points  respectively,  (c)  and  (d)  balance 
of  the  incoming  and  outgoing  heat  fluxes  for  the  grid  of  4 1  X  101  points 
and  81  X  181  points  respectively. 


required  to  maintain  mass  conservation  in  the  molten  zone  decreases  from  the  normal  gravity 
case  shown  in  Figure  5.11.  Figure  5.16  compares  the  temperature  profiles  between  the  41 
X  101  and  the  81  X  181  grids  showing  grid  independence  of  the  solution.  It  is  also  noted, 
from  comparing  Figures  5.9(c  -  d)  and  5. 13(c  -  d)  that  the  microgravity  situation  causes  a 
larger  heat  flux  to  be  transported  through  the  zone.  This  is  due  to  the  larger  meniscus  surface 
area  in  the  vicinity  of  the  heater. 
5.6.1.2  Buoyancy  Driven  Convection 

Figures  5.17  -  5.20  show  the  effect  of  buoyancy  driven  convection  on  the  solution 
characteristics.  Thermocapillary  convection  has  been  suppressed  by  imposing  a  zero  shear 
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Hm  =  0.327  cm, 
(i)   Hc  =  0.234  cm. 


Hm  =  0.330  cm, 
(ii)  Hc  =  0.239  cm. 


Figure  5.14  Computed  grid  and  boundary  shape  for  the  conduction  solution  at  Og.  (i) 
Grid  of  41  X  101  points,  (ii)  Refined  grid  of  81  X  181  points.  Only  the  phase  change 
region  and  its  immediate  vicinity  is  shown  since  the  actual  height  of  the  rod  is  40  times 
the  radius.  The  zone  heights  agree  to  1%  at  the  surface  and  2%  at  the  centerline. 


Figure  5.15  Isotherms  for  the  conduction  solution  at  Og.  (i)  Grid  of  41  x  101  points 
and  (ii)  Refined  grid  of  81  X  181  points.  The  Tmax  agree  to  0.1%. 
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Temperature  profiles  for  conduction  solution  and  deforming  meniscus  at  Og 
Surface  Temperature  Profiles  Centerline  Temperature  Profiles 

8  i  '  1  ■  1  '  1  1  8 


0 


8 


f 
5 
= 


T  :  81  x  181 
T. :  41 x 101 


300       400       500  600 

(a) 

41  x  101  grid 


0 

300       400  500 

(c)  T(K) 


600 


300 

(b) 


300 

(d) 


T  :  81 x 181 
T. :  41  x  101 


400  500 


81  x  181  grid 


600 


400  500 

T(K) 


600 


Figure  5.16  (a)  and  (b)  Temperature  profiles  showing  comparison 
between  coarse  grid  of  41  X  101  points  and  fine  grid  of  81  X  181 
points  along  the  entire  height  of  the  domain,  (a)  Temperature  profile 
along  the  surface,  (b)  Temperature  profile  along  the  centerline. 
(c)  and  (d)  Comparison  of  Temperatures  along  the  surface  and  the 
centerline.  (c)  41  X  101  grid,  (d)  81  X  181  grid. 


stress  condition  on  the  free  surface.  Convergence  was  declared  when  the  residuals  were  less 
than  10-5  for  the  momentum  equations  and  10-7  for  the  energy  equation.  The  1-norm  of  the 
solution  was  less  than  10-5  for  all  the  dependent  variables.  Figure  5.17  (a)  and  (b)  show  the 
convergence  history  of  the  zone  lengths -Hm  and  Hc- on  the  41  X  101  and  81  X  181  grids, 
respectively.  Figure  5.17(c)  and  (d)  show  the  history  of  the  heat  flux  balance  in  the  domain 
on  the  41  X  101  and  81  X  181  grid  systems,  respectively.  The  heat  flux  balance  on  the  41 
X  101  grid  is  satisfied  to  within  0.2%  but  on  the  81  X  181  grid  there  is  a  2%  imbalance. 
However,  this  is  adequate  for  the  purpose  of  accuracy  assessment  and  comparison  with  the 
solutions  of  Lan  (1994).  The  streamfunction  is  shown  in  Figure  5.18.  The  maximum 
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Solution  history  for  natural  convection  and  deforming  meniscus 
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Figure  5.17  Convergence  history  for  the  natural  convection  solution,  (a) 
and  (b)  height  of  the  zone  at  the  surface  and  at  the  centerline  for  the  grid  of 
41  x  101  points  and  81  X  181  points  respectively,  (c)  and  (d)  balance  of 
the  incoming  and  outgoing  heat  fluxes  for  the  grid  of  41  X  101  points 
(0.2%)  and  81  X  181  points  (2%)  respectively. 


streamfunction,  Wmax,  differs  by  less  than  2%.  The  present  solution  on  the  41  X  101  grid 
differs  by  8%  from  the  maximum  streamfunction  reported  in  Lan  (1994).  Figure  5. 19  shows 
the  corresponding  isotherms.  The  maximum  temperature,  Tmax,  occuring  on  the  free  surface 
differs  by  0.05%  between  the  solutions  and  0.4%  from  the  solution  of  Lan  (1994).  Figure 
5.20  (a)  and  (b)  show  identical  temperature  profiles  for  the  solutions  on  41  X  101  and  81 
X  181  grids  indicating  grid  independence  of  the  solution.  Figures  5.20(c)  and  (d)  may  be 
contrasted  with  the  conduction  case  illustrated  in  Figures  5.12(c)  and  (d)  showing  the 
decrease  in  Tmax  due  to  the  effectiveness  of  convection  in  heat  transport.  The  effect  of  natural 
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(i) 

*Pmax  =  2.706  X  10"3g/s, 
AW  =  1.87  X  10"4  g/s 


(ii) 


Wmax  =  2.654  X  10-Jg/s 
A*P=  1.84  X  l(Hg/s 


Figure  5.18.  Streamfunction  and  phase  interfaces  on  (i)  41  X  101  and  (ii)  81  X 
181  grids,  (i)  and  (ii)  differ  by  less  than  2%.  The  arrow  indicates  the  direction  of 
the  convection  at  the  free  surface.  The  41  X  101  grid  solution  differs  from  Lan 
(1994)  by  8%. 


(i) 

AT  =  3.99  K,  Tmax  =  601.03  K, 


(ii) 

AT  =  4.07  K,  Tmax  =  600.76  K, 


min 


=  301.23  K 


min 


=  300.93  K 


Hm  =  0.344  cm,  Hc  =  0.257  cm 
PL  =  995.46  dynes/cm2 


Hm  =  0.335  cm,  Hc  =  0.244  cm 
PL  =  995.31  dynes/cm2. 


Figure  5.19  Isotherms  on  (i)  41  x  101  and  (ii)  81  X  181  grids.  Tmax  differs  by 
0.05%.  The  41  X  101  grid  solution  differs  from  the  solution  of  Lan  (1994)  by  0.4%. 


149 


Temperature  profiles  for  natural  convection  and  deforming  meniscus 
Surface  Temperature  Profiles  Centerline  Temperature  Profiles 
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Figure  5.20  (a)  and  (b)  Temperature  profiles  showing  comparison  between  coarse 
grid  of  41  X  101  points  and  fine  grid  of  81  X  181  points  along  the  entire  height  of 
the  domain,  (a)  Temperature  profile  along  the  surface,  (b)  Temperature  profile  along 
the  centerline. 


convection  is  to  cause  the  heated  melt  to  rise  along  the  free  surface,  transfer  heat  to  the 
melt-feed  interface,  and  sink  at  the  centerline  towards  the  melt-crystal  interface.  This 
causes  the  melt-feed  interface  to  acquire  a  more  distorted  shape  compared  with  the 
conduction  solution,  and  the  melt-crystal  interface  to  become  flatter  over  most  of  the  region 
except  close  to  the  meniscus.  The  distortion  of  the  isotherms  and  the  upward  shift  in  the 
location  of  Tmax  is  due  to  the  convective  action.  Since  the  Prandtl  number  is  high,  even  a 
modest  Grashof  number  results  in  a  significant  impact  on  the  temperature  distribution.  The 
higher  temperature  gradients  evident  at  the  base  of  the  zone  also  indicate  that  the  heating 
takes  place  close  to  the  bottom  of  the  zone.  The  zone  is  translated  upward  as  compared  with 
the  conduction  solution  due  to  the  enhancement  of  heat  transfer  through  the  upper  half  of 
the  domain  and  its  suppression  throught  the  lower  half. 

5.7  Thermocapillary  Convection 
In  the  following  set  of  calculations,  the  effect  of  thermocapillary  convection  is  taken 
into  account.  As  will  be  demonstrated,  the  Marangoni  effect  has  a  significant  impact  on  the 


150 


temperature  distribution  and  the  phase  interfaces  within  the  melt  zone.  The  actual  Marangoni 
numbers  in  this  particular  system  are  in  the  region  of  5000,  however,  numerical  stability 
constrained  the  calculations  to  a  Marangoni  number  of 500.  Grid  independence  and  accuracy 
assessments  have  been  carried  out.  Computations  have  been  carried  out  to  assess  the  impact 
of  meniscus  deformation  on  the  solid  -  liquid  interface  shape  and  zone  size.  Calculations 
have  also  been  conducted  to  simulate  float  zone  growth  under  zero-gravity  conditions  to 
simulate  float  zone  processing  in  orbit.  In  all  subsequent  calculations,  the  convergence 
criteria  has  been  set  to  less  than  10-4  for  the  momentum  equations  and  to  less  than  10-6  for 
the  energy  equations  while  ensuring  that  the  1-norm  of  the  solution  was  less  than  10-5.  These 
criteria  either  meet  or  exceed  the  convergence  criteria  specified  in  the  literature  cited 
previously. 

5.7. 1    Normal  Gravity 
5.7.1.1  Cylindrical  meniscus 

A  conduction  calculation  with  the  imposed  cylindrical  meniscus  was  used  as  the  initial 
guess.  Figure  5.21  shows  the  convergence  history  for  calculations  on  the  41  X  101  grid  as 
well  as  the  81  X  181  grid.  Figure  5.21(a)  and  (b)  show  the  history  of  the  zone  height  and 
Figure  5.21(c)  and  (d)  show  the  heat  flux  balance  within  the  domain.  Figure  5.22(a)  shows 
the  streamfunction  and  Figure  5.22(b)  shows  the  computed  isotherms.  The  streamfunction 
pattern  shows  the  interaction  between  natural  and  thermocapillary  convection.  Two  large 
convection  rolls  are  apparent  in  the  upper  half  of  the  domain  where  the  thermocapillary 
convection  aids  natural  convection.  In  the  lower  part  of  the  zone,  two  small  counter-rotating 
loops  are  formed  where  the  thermocapillary  effect  opposes  natural  convection.  The  isotherm 
plots  show  that  the  location  of  Tmax  is  displaced  upwards  and  its  value  is  less  as  compared 
with  the  conduction  case.  The  distortion  of  the  phase  interfaces  is  clearly  evident  especially 
away  from  the  centerline.  Figure  5.23(a)  and  (b)  provide  a  comparison  of  the  temperature 
profiles  on  the  two  grid  resolutions  considered.  The  temperature  distribution  is  almost 
identical  on  the  two  grids. 
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Solution  history  for  Ma  =  500,  lg,  and  cylindrical  meniscus 
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Figure  5.21  Convergence  history  for  the  thermocapillary  convection 
solution  Ma  =  500,  at  lg  and  with  a  fixed  cylindrical  meniscus,  (a)  and  (b) 
height  of  the  zone  at  the  surface  and  at  the  centerline  for  the  grid  of  41  x 
101  points  and  81  X  181  points  respectively,  (c)  and  (d)  balance  of  the 
incoming  and  outgoing  heat  fluxes  for  the  grid  of  41  X  101  points  (1.7%) 
and  81  X  181  points  (7%)  respectively. 


5.7.1.2  Deformed  meniscus 

The  previous  calculation  is  used  as  the  starting  guess  and  now  the  meniscus  is  allowed 
to  deform  according  to  the  requirement  of  mass  conservation  and  the  Young-Laplace 
equation.  Figure  5.24  shows  the  convergence  history  of  the  solution.  Figure  5.24(a)  and  (b) 
show  the  convergence  history  of  the  zone  height  and  an  initially  oscillatory  convergence  path 
can  be  observed  which  damps  out  as  the  iterations  continue.  Figure  5.24(c)  and  (d)  show  the 
heat  flux  balance  for  the  solutions  on  the  41  X  101  grid  and  the  8 1  X  181  grid  respectively. 
For  the  41  X  101  grid,  a  heat  flux  imbalance  of  4%  is  observed  which  reduces  to  0.8%  with 
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(i) 


Wmax=  1.733  x  l(Hg/s,(+) 
A*P  =  1.20  X  1(H  g/s 
^1  = -4.433  X  10"4  g/s 


(ii) 


^max 
A*P=  1.55 
=  -4.840 


=  2.242  X  10-3  g/s 
X  10^  g/s 
X  1(H  g/s 


(a)  Streamfunction  and  phase  interfaces  on  (i)  41  X  101  and  (ii)  81  X  181  grids. 
Arrows  show  the  direction  of  convection  cells. 


(i) 


(ii) 


AT  =  3.90  K,  Tmax  =  596.97  K, 

Tmin  =  300.85  K 

Hm  =  0.331  cm,  Hc  =  0.166  cm 


AT  =  4.04  K,  Tmax  =  597.81  K, 

Tmin  =  300.94  K 

Hm  =  0.365  cm,  Hc  =  0.223  cm 


(b)  Isotherms  on  (i)  41  x  101  and  (ii)  81  X  181  grids. 


Figure  5.22  Solution  for  Ma  =  500,  lg,  and  fixed  cylindrical  meniscus. 
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Temperature  profiles  for  Ma  =  500,  lg,  and  cylindrical  meniscus 
Surface  Temperature  Profiles  Centerline  Temperature  Profiles 


Figure  5.23  (a)  and  (b)  Temperature  profiles  showing  comparison  between 
coarse  grid  of  41  X  101  points  and  fine  grid  of  8 1  X  181  points  along  the  entire 
height  of  the  domain,  (a)  Temperature  profile  along  the  surface,  (b)  Temperature 
profile  along  the  centerline. 


grid  refinement.  Figure  5.25(a)  shows  the  streamfunction  on  the  41  X  101  grid  and  the  81 
X  181  grid.  Figure  5.25(b)  shows  the  corresponding  isotherms.  Allowing  the  meniscus  to 
deform  increases  the  free  surface  area  thereby  increasing  the  convection  strength  of  all  the 
four  convection  cells  in  the  domain.  Since  the  maximum  Temperature,  Tmax,  is  virtually 
unaffected,  this  implies  a  larger  heat  flux  through  the  domain  resulting  in  a  taller  melt  zone. 
The  increased  strength  of  the  upper  convection  rolls  also  causes  the  melt-crystal  interface 
at  the  lower  end  of  the  melt  zone  to  become  flatter  at  the  centerline,  but  acquire  increased 
distortion  close  to  the  free  surface.  The  maximum  streamfunction  differs  by  9.7%  between 
the  solutions  on  the  grid  of  41  X  101  points  and  the  grid  of  81  X  181  points.  Figure  5.26(a) 
and  (b)  compare  the  temperature  profiles  of  the  two  solutions  which  are  indistinguishable. 
Figure  5.26(c)  and  (d)  show  the  upward  displacement  of  the  location  of  Tmax  as  compared 
with  the  conduction  solutions  and  the  establishment  of  the  one  dimensional  temperature 
distribution  in  the  feed  and  crystal  regions  away  from  the  melt  zone. 
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Solution  history  for  Ma  =  500,  lg,  and  deforming  meniscus 
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Figure  5.24  Convergence  history  for  the  thermocapillary  convection 
solution  Ma  =  500,  at  lg  and  with  computed  meniscus  shape,  (a)  and  (b) 
height  of  the  zone  at  the  surface  and  at  the  centerline  for  the  grid  of  41  X 
101  points  and  81  X  181  points  respectively,  (c)  and  (d)  balance  of  the 
incoming  and  outgoing  heat  fluxes  for  the  grid  of  41  X  101  points  (4%) 
and  81  X  181  points  (0.8%)  respectively. 


5.7.1.3  Zero  Gravity 

5.7.1.4  Cylindrical  meniscus 

A  conduction  solution  with  an  imposed  cylindrical  meniscus  was  used  as  the  initial 
guess.  Figure  5.27  shows  the  convergence  history  of  the  solution.  Figure  5.27(a)  and  (b) 
show  the  history  of  the  zone  height  for  the  solutions  on  the  4 1  x  101  grid  and  the  8 1  X  181 
grid  respectively.  Figure  5.27(c)  and  (d)  show  the  corresponding  heat  flux  balance.  Figure 
5.28(a)  shows  the  streamfunction  and  Figure  5.28(b)  the  isotherm  distribution.  It  is 
immediately  evident  that  the  solution  is  symmetrical  about  y  =  4  cm  which  is  the  heater 
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(i) 


(ii) 


*I>max=  2.880  X  10~3g/s, 
A^=  1.99  X  10"4  g/s 
Wl  =  ^.975  x  1(H  g/s 

(a)  Streamf unction  and  phase  interfaces  on  (i)  41  x  101  and  (ii)  81  X  181  grids. 


Wmax  =  2.623  X  10"3g/s 
A^=  1.81  X  lO^g/s 
^1  = -5.160  x  10-*  g/s 


(i) 


(ii) 


AT  =  4.01  K,  Tmax  =  598.66  K, 
Train  =  301.33  K 
Hm  =  0.394  cm,  Hc  =  0.260  cm. 
PL=  1011.47  dynes/cm2. 


AT  =  4.07  K,  Tmax  =  597.87  K, 


min 


=  300.94K 


Hm  =  0.372  cm,  Hc  =  0.227  cm 
PL=  1011.52  dynes/cm2. 


(b)  Isotherms  on  (i)  41  X  101  and  (ii)  81  X  181  grids. 


Figure  5.25  Solution  for  Ma  =  500,  lg,  and  computed  meniscus. 
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Temperature  profiles  for  Ma  =  500,  lg,  and  deforming  meniscus 
Surface  Temperature  Profiles  Centerline  Temperature  Profiles 


Figure  5.26  (a)  and  (b)  Temperature  profiles  showing  comparison  between  coarse 
grid  of  41  X  101  points  and  fine  grid  of  81  X  181  points  along  the  entire  height  of 
the  domain,  (a)  Temperature  profile  along  the  surface,  (b)  Temperature  profile  along 
the  centerline. 


location.  Four  counter-rotating  flow  loops  are  evident  in  the  melt  zone  and  these  convection 
rolls  cause  the  isotherm  distribution  to  be  pulled  towards  the  free  surface.  Convection 
suppresses  heat  transfer  to  the  centerline  causing  the  characteristic  convexity  of  the  phase 
interfaces  towards  the  melt.  The  maximum  streamfunction  is  lower  than  that  in  the  normal 
gravity  case,  but  the  flattening  effect  of  natural  convection  on  the  melt-crystal  interface  is 
absent.  Figure  5.29(a)  and  (b)  compare  the  temperature  profiles  on  the  two  grids  and  it  may 
be  observed  that  the  two  solutions  are  nearly  identical. 
5.7.1.5  Deformed  meniscus 

The  previously  computed  solution  was  used  as  the  initial  guess  and  the  meniscus  was 
allowed  to  deform  in  response  to  the  needs  of  global  mass  conservation  and  the 
Young-Laplace  equation.  For  this  case,  momentum  residuals  of  less  than  5  X  10-4  could 
not  be  achieved  on  the  grid  of  41  X  101  points.  Figure  5.30  shows  the  convegence  history 
of  the  solution.  Figure  5.31(a)  shows  the  streamfunction  and  Figure  5.31(b)  shows  the 
isotherm  distribution  on  the  grid  of  41  x  101  and  81  X  181  points.  It  is  observed  that  the 
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Solution  history  for  Ma  =  500,  Og,  and  fixed  cylindrical  meniscus 
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Figure  5.27.  Convergence  history  for  the  thermocapillary  convection 
solution  Ma  =  500,  at  Og  and  with  fixed  cylindrical  meniscus,  (a)  and  (b) 
height  of  the  zone  at  the  surface  and  at  the  centerline  for  the  grid  of  41  X 
101  points  and  81  x  181  points  respectively,  (c)  and  (d)  balance  of  the 
incoming  and  outgoing  heat  fluxes  for  the  grid  of  41  X  101  points  (8.6%) 
and  81  X  181  points  (4.6%)  respectively. 


convection  rolls  increase  in  strength  since  the  meniscus  deformation  results  in  a  larger 
surface  area  over  which  the  surface  tension  gradients  can  act.  These  factors  also  result  in  a 
larger  heat  flux  being  transported  through  the  domain  and  the  height  of  the  melt  zone  to 
increase.  The  maximum  streamfunction  differs  by  6.6%  between  the  solution  on  the  41  X 
101  grid  and  81  X  181  grid.  Figure  5.32(a)  and  (b)  show  that  the  temperature  profiles 
obtained  on  the  two  grids  are  identical. 
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(i) 

AT  =  3.65  K,  Tmax  =  600.18  K, 

Tmin  =  298.00  K 

Hm  =  0.374  cm,  Hc  =  0.199  cm 


(") 

AT  =  4.14  K,  Tmax  =  599.50  K, 


min 


=  298.00  K 


Hm  =  0.356  cm,  Hc  =  0. 173  cm 


(b)  Isotherms  on  (i)  41  X  101  and  (ii)  81  X  181  grids. 


Figure  5.28  Solution  for  Ma  =  500,  Og,  and  fixed  cylindrical  meniscus. 
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Temperature  profiles  for  Ma  =  500,  Og,  and  fixed  cylindrical  meniscus 
Surface  Temperature  Profiles  Centerline  Temperature  Profiles 


Figure  5.29  (a)  and  (b)  Temperature  profiles  showing  comparison  between  coarse 
grid  of  41  X  101  points  and  fine  grid  of  81  X  181  points  along  the  entire  height  of 
the  domain,  (a)  Temperature  profile  along  the  surface,  (b)  Temperature  profile  along 
the  centerline. 


5.8  Conclusion 

A  moving  boundary  technique  has  been  developed  in  the  context  of  the  SIMPLE 
algorithm  and  the  enthalpy  method  for  the  computation  of  convective  heat  transfer  and 
solidification  in  the  float  zone.  A  general  multi-block  technique  for  algebraic  grid 
generation  has  been  developed  to  fit  curved  free  surfaces  and  adaptively  enclose  the  phase 
change  regions  in  solidification  processes.  Underrelaxation  procedures  have  been  developed 
to  ensure  the  stability  of  the  numerical  procedure.  Care  has  been  taken  to  ensure  that  grid 
quality  is  maintained  after  the  grid  generation  procedure  with  appropriate  error  recovery  in 
case  of  failure  of  the  grid  generation  procedure. 

Computations  have  been  carried  out  to  assess  the  robustness,  stability  and  accuracy  of 
the  convective  heat  transfer  and  phase  change  computations  via  grid  refinement  and 
comparison  with  the  solutions  of  Lan  (1994).  The  solutions  of  Lan  (1994)  were  based  on 
front  tracking  rather  than  the  enthalpy  method,  used  the  streamfunction-vorticity 
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Solution  history  for  Ma  =  500,  Og,  and  deforming  meniscus 
41  x  101  grid  81  x  181  grid 
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Figure  5.30  Convergence  history  for  the  thermocapillary  convection 
solution  Ma  =  500,  at  Og  and  with  computed  meniscus  shape,  (a)  and 
(b)  height  of  the  zone  at  the  surface  and  at  the  centerline  for  the  grid 
of  41  X  101  points  and  81  X  181  points  respectively,  (c)  and  (d) 
balance  of  the  incoming  and  outgoing  heat  fluxes  for  the  grid  of  41  X 
101  points  (6%)  and  81  X  181  points  (2.7%)  respectively. 


formulation  of  the  Navier-Stokes  equations  as  opposed  to  the  primitive  variables 
formulation  used  in  the  present  work,  and  employed  first  order  upwind  schemes  for  the 
convection  schemes  rather  than  the  second  order  central  difference  scheme  employed  in  this 
study.  It  is  found  that  for  the  conduction-capillarity  and  buoyancy  driven  convection 
calculations,  grid  independent  solutions  can  be  obtained  even  on  41  X  101  grids.  However, 
when  thermocapillary  convection  is  taken  into  account,  grid  refinement  from  the  41  X  101 
grid  to  the  81  X  181  grid  exhibits  quantitative  differences  although  not  qualitative 
differences.  Good  comparison  has  been  obtained  with  the  solutions  of  Lan  (1994)  for  the 
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(i) 


(ii) 


Wmax=  1.097  X  10-3g/s 


Wmax=  1.170  X  10-3  g/s, 

AW  =  8.07  X  10-5  g/s  AW  =  7.57  X  10~5  g/s 

(a)  Streamfunction  and  phase  interfaces  on  (i)  41  X  101  and  (ii)  81  X  181  grids. 


(i) 


AT  =  3.63  K,  Tmax  =  600.89  K, 


min 


=  298.00  K 


Hm  =  0.406  cm,  Hc  =  0.219  cm 
PL  =  642.25  dynes/cm2 


(ii) 


AT  =  4.23  K,  Tmax  =  599.95  K, 


min 


=  300.9  IK 


Hm  =  0.381  cm,  Hc  =  0.183  cm 
Pl  =  654.94  dynes/cm2. 


(b)  Isotherms  on  (i)  41  X  101  and  (ii)  81  X  181  grids. 


Figure  5.31  Solution  for  Ma  =  500,  Og,  and  computed  meniscus. 
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Temperature  profiles  for  Ma  =  500,  Og,  and  deforming  meniscus 
Surface  Temperature  Profiles  Centerline  Temperature  Profiles 
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Figure  5.32  (a)  and  (b)  Temperature  profiles  showing  comparison  between 
coarse  grid  of  41  X  101  points  and  fine  grid  of  81  X  181  points  along  the  entire 
height  of  the  domain,  (a)  Temperature  profile  along  the  surface,  (b)  Temperature 
profile  along  the  centerline. 


conduction-capillarity  and  natural  convection  solutions  with  similar  boundary  shapes  and 
grid  distributions. 

Computations  have  been  carried  out  to  assess  the  impact  of  free  surface  deformation  as 
well  as  microgravity  environment  on  the  float  zone.  The  free  surface  deformation  results  in 
a  larger  surface  area  at  the  meniscus  resulting  in  a  higher  strength  of  thermocapillary 
convection  as  well  as  an  increase  in  the  heat  flux  through  the  system.  Thus,  for  the  parameters 
and  regimes  considered  in  this  study,  the  free  surface  deformation  has  a  quantitative  but  not 
qualitative  impact  on  the  melt  zone.  The  microgravity  environment,  on  the  other  hand, 
results  in  a  suppression  of  buoyancy  induced  convection  as  well  as  a  significant  impact  on 
the  meniscus  shape.  However,  at  the  operating  parameters  considered  in  this  study,  the 
microgravity  environment  does  not  offer  any  significantly  reduced  convection  since  the 
surface  tension  gradient  is  always  present  due  to  the  thermal  inhomogeneties  along  the 
meniscus. 
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The  difference  between  the  present  study  and  the  previous  studies  lies  in  the  use  of  the 
enthalpy  method  as  opposed  to  interface  tracking  and  the  use  of  a  flexible,  multi-block  grid 
generation  procedure  that,  in  conjunction  with  the  use  of  the  enthalpy  method,  can  adapt  to 
a  wide  range  of  zone  lengths  including  the  fully  solid  case.  The  procedure  has  been  designed 
to  accommodate  unsteady  computations,  to  be  presented  in  a  subsequent  study,  which  can 
use  a  fully  solid  rod  as  the  initial  condition.  A  front  tracking  procedure  as  employed  by  Lan 
( 1 994)  or  Lan  and  Kou  ( 1 99 1 )  cannot  handle  the  small  zone  sizes  that  will  be  present  initially. 
The  use  of  the  enthalpy  method  in  this  study  makes  it  unnecessary  to  fit  the  solid-liquid 
interfaces  giving  additional  flexibility.  This  study  also  uses  higher  order  convection  schemes 
-  the  second  order  central  difference  scheme  as  opposed  to  the  first  order  upwind  scheme 
used  in  the  previous  studies.  Also,  grid  refinement  studies  have  been  carried  out  employing 
substantially  higher  grid  resolutions  than  the  previous  studies.  The  study  indicates  that 
higher  grid  resolutions,  such  as  those  considered  in  this  study,  are  needed  as  the  strength  of 
thermocapillary  convection  increases  in  the  present  float  zone  system. 


CHAPTER  6 
CONCLUSION 


6. 1  Summary  and  Conclusions 
A  computational  tool  has  been  developed  to  simulate  fluid  flow,  heat  transfer,  moving 
free  surfaces,  and  phase  change  dynamics  in  the  context  of  problems  involving  solidification 
and  single  crystal  growth. 

( 1 )  The  SIMPLE  algorithm  in  body  fitted  coordinates  has  been  used  as  a  starting  point  for 
algorithm  development. 

(2)  A  previously  developed  free  energy  minimization  model  has  been  used  to  compute 
meniscus  shapes  based  on  the  static  Young-Laplace  equation  for  problem  of  float  zone 
growth  of  NiAl.  The  free  energy  criterion  has  been  verified  with  available 
experimental  data  and  good  agreement  obtained.  The  free  energy  criterion  has  then 
been  used  to  generate  meniscus  shapes  for  the  float  zone  problem.  When  dynamical 
effects  of  melt  convection  are  taken  into  account,  it  is  found  that  the  Capillary  numbers 
are  very  small  compared  to  the  Bond  number  and  the  viscous  normal  stress  is  too  small 
to  have  any  significant  impact  on  the  meniscus  shape  thus  justifying  the  static 
meniscus  shape  used  in  the  solidification  calculations. 

(3)  A  moving  boundary  computational  procedure  has  been  formulated  in  order  to  couple 
the  free  surface  computation  with  the  solid/liquid  interface  dynamics  and  melt 
convection.  An  algebraic  grid  generation  procedure  using  transfinite  interpolation  has 
been  used  to  regenerate  the  grid  at  every  iteration  following  the  update  of  the  meniscus 
shape  from  the  normal  stress  balance  implied  by  the  dynamic  Young-Laplace 
equation. 
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(4)  A  variety  of  accurate  convection  schemes,  the  second  order  central  difference  scheme, 
the  second  order  upwind  scheme,  and  a  'mixed'  formulation  consisting  of  the  second 
order  upwind  scheme  for  the  momentum  equations  and  the  lower  order  hybrid  scheme 
for  the  energy  equation  have  been  tested  against  the  benchmark  problem  of  natural 
convection  of  air  in  a  square  cavity.  It  has  been  found  that  the  'mixed'  formulation 
suffered  no  degradation  of  accuracy  in  spite  of  the  formal  lower  order  of  accuracy  even 
for  high  Rayleigh  numbers. 

(5)  Time  dependent  simulations,  using  the  second  order  central  difference  scheme  for  the 
convection  terms  and  the  first  order  backward  Euler  method  for  the  transient  terms, 
have  been  conducted  for  the  problem  of  transient  natural  convection  around  an 
enclosed  vertical  channel,  with  application  to  electronics  cooling.  Care  has  been  taken 
to  maintain  good  accuracy  in  space  and  time  by  means  of  grid  and  time  step 
refinement. 

(6)  Steady,  thermocapillary  convection  in  an  idealized  model  for  float  zone  growth  of 
NiAl  single  crystals  has  been  simulated.  In  the  simplified  model,  the  liquid  bridge  has 
been  considered  to  be  suspended  between  two  fixed  flat  plates  and  the  meniscus  shape 
satisfies  the  static  Young-Laplace  equation.  Consequently,  the  meniscus  shape  is 
imposed  in  this  static  model  and  is  not  computed  as  part  of  the  solution.  It  is  found  that 
the  isotherm  distribution  for  high  Marangoni  number  convection  differs  considerably 
from  the  heat  conduction  solution  obtained  by  neglecting  convective  effects 
altogether.  It  is  also  determined  that  thermocapillary  convection  dominates  over 
natural  convection  even  at  1-g  conditions. 

(7)  A  single  region  enthalpy  formulation  has  been  developed  for  the  numerical 
computation  of  phase  change  processes  on  a  fixed  grid.  Two  alternative 
implementations  of  the  enthalpy  method  have  been  developed  -  the  T-based  method 
and  the  H-based  method.  The  one  dimensional  Stefan  problem  has  been  used  to  test 
the  accuracy  and  efficiency  of  the  two  methods.  It  has  been  found  that  while  both 
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methods  are  equally  accurate,  the  T-based  method  is  faster  and  computationally 
efficient  over  a  range  of  Stefan  numbers  and  mushy  zone  thickness.  The  H-based 
method  is  more  robust  over  a  wide  range  of  parameters  and  works  even  when  the  size 
of  the  mushy  zone  shrinks  to  zero. 

(8)  The  enthalpy  porosity  model  has  also  been  tested  on  a  two  dimensional  melting 
problem  where  the  heat  transfer  is  dominated  by  natural  convection  in  the  melt. 
Qualitative  agreement  has  been  obtained  with  available  experimental  data. 

(9)  The  simplified  float  zone  model  developed  for  single  crystal  growth  of  Ni  Al  has  been 
improved  by  relaxing  the  restriction  of  flat  solid/liquid  interfaces  and  employing  a 
moving  grid  procedure  to  track  the  deforming  free  surface  of  the  float  zone.  The 
enthalpy  porosity  method  has  been  used  to  couple  the  phase  change  computations  with 
thermocapillary  and  natural  convection.  Computations  have  been  carried  out  to  assess 
the  interaction  between  thermocapillary  convection  and  buoyancy  induced  convection 
and  the  impact  of  free  surface  deformation  by  considering  normal  gravity  as  well  as 
microgravity  conditions. 

6.2  Future  Work 

6.2.1    Node  Jumping 

The  enthalpy  method  is  prone  to  the  node  jumping  problem  even  for  the  unsteady  one 
dimensional  conduction  problem  as  is  observed  in  Fig.  4.5  -  4.7.  Different  researchers  have 
attempted  to  deal  with  this  problem  by  diagnosing  that  the  node  jumping  is  due  to  the  nodal 
temperature  of  the  phase  change  control  volume  being  held  constant  at  Tm  even  when  the 
interface  has  moved  away  from  the  nodal  point.  Voller  and  Cross  (1981)  iteratively  adjust 
the  time  step  so  that  the  interface  is  always  located  at  the  nodal  points.  Date  (1992)  has 
attempted  to  correct  the  nodal  temperatures  by  estimating  the  location  of  the  interface  in  a 
one  dimensional  context  to  obtain  a  smooth  motion  of  the  interface.  However,  such 
corrections  are  not  straightforward  to  obtain  in  a  multidimensional  problem  and  can  be 
implemented  in  many  different  ways.  In  the  current  study  of  the  enthalpy  method  shown  in 
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Fig.  4.5  -  4.7  it  is  observed  that  the  node  jumping  problem  decreases  as  the  thickness  of  the 
mushy  zone  increases;  nevertheless,  it  is  worthwhile  to  cure  the  node  jumping  problem  in 
a  multidimensional  context  to  obtain  better  time  accuracy. 

6.2.2  Multi-Block  Approach 

A  multi-block  approach  (Wright  and  Shyy  1992,  Shyy  et.  al.  1994)  to  the  float  zone 
problem  would  result  in  better  accuracy  and  computational  economy.  For  instance,  the  single 
region  formulation  implies  that  all  the  field  equations  are  solved  in  both  the  phases  which 
is  wasteful  of  computer  time  and  storage.  In  a  multi-block  approach,  only  the  heat 
conduction  equation  needs  to  be  solved  in  the  crystal  and  the  feed  regions;  also,  fewer  grid 
points  in  the  radial  direction  can  be  used  in  these  regions  since  the  temperature  profiles  vary 
slowly  and  are  predominantly  one  dimensional  in  these  regions.  This,  coupled  with  grid 
adaptation  towards  the  mushy  zone  and  regions  of  high  velocity  gradients  (Shyy  1987),  will 
allow  more  grid  points  to  be  located  in  the  phase  change  and  moving  boundary  region, 
resulting  better  accuracy  and  economy. 

6.2.3  Feed/Crystal  Rotation 

Some  researchers  (Lan  and  Kou  1991c,  1992b,  1993b,  Muhlbauer  et.  al.  1995)  have 
determined  that  rotation/counterrotation  of  the  feed  and/or  the  crystal  can  result  improved 
heat  transfer  to  the  centerline  of  an  otherwise  axisymmetric  melt  due  to  the  swirl  induced 
in  the  melt  resulting  in  shorter  zones  and  enabling  larger  diameter  crystals  to  be  grown. 

6.2.4  Dopant  Transport 

Dopant  transport  and  concentration  play  an  important  role  in  determining  the 
solid/liquid  interface  shape  and  stability.  Inclusion  of  solute  transport  and  segregation 
models  will  make  it  more  complete  and  useful. 

6.2.5  Electromagnetic  Heating 

A  more  complete  model  of  the  float  zone  process  can  be  constructed  by  including 
electromagnetic  effects.  For  electrically  conducting  materials,  induction  heating  is  a 
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common  means  of  heating  the  zone.  Associated  electromagnetic  levitation  forces  also  help 
to  stabilize  the  meniscus  against  gravity  enabling  bigger  zone  sizes  under  1-g  conditions. 
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